Seismic characteristics of the 15 February 2013 bolide explosion in Chelyabinsk, Russia

The seismological characteristics of the 15 February 2013 Chelyabinsk bolide explosion are investigated based on seismograms recorded at 50 stations with epicentral distances ranging from 229 to 4324 km. By using 8–25 s vertical‐component Rayleigh waveforms, we obtain a surface‐wave magnitude of 4.17±0.31 for this event. According to the relationship among the Rayleigh‐wave magnitude, burst height and explosive yield, the explosion yield is estimated to be 686 kt. Using a single‐force source to fit the observed Rayleigh waveforms, we obtain a single force of 1.03×1012 N, which is equivalent to the impact from the shock wave generated by the bolide explosion.

Earth's atmosphere is frequently impacted by asteroids that can result in bolides. Over 500 separate events have occurred in the last 20 years (http://www.jpl.nasa.gov/news/news.php?release= 2014-397). Most of the asteroids disintegrate in the atmosphere, due to their small sizes, without doing any harm to the Earth's fauna and flora. But some big events, such as the 1908 Tunguska event in Russia, have been devastating. Large bolide explosions can create shock waves that cause considerable damage upon impact with Earth's surface, and even trigger earthquakes. Seismic records of the 15 February 2013 Chelyabinsk bolide explosion are investigated quantitatively in order to improve estimations of the threat to Earth posed by impact of large asteroids. Antolik et al. (2014) calculated a Rayleigh-wave magnitude of 4.06±0.22 and, by assuming a coupling coefficient of 0.001% between the yield and , obtained an estimate of approximately 300 kt for the energy released from the bolide blast.
Furthermore, Heimann et al. (2013) investigated the seismic source parameters of the Chelyabinsk event through full waveform fitting based on an isotropic atmospheric airburst model, and Antolik et al. (2014) obtained the source mechanism parameters for the explosion based on a full moment tensor solution.
When shock waves propagate through the air, their energy decays rapidly with distance, and the manner in which air waves impact on the ground can be represented as distributed forces. Considering that the area influenced by the shock waves of this event was limited, the impact can be simplified as a single force. Therefore, in this study, we use a single-force model to investigate the 2013 bolide explosion.

Data
We collected seismograms recorded at 50 Global Seismic Network (GSN) stations located at regional and teleseismic distances between 229 km and 4320 km from the epicenter of the bolide explosion. All stations are equipped with broadband instruments with nearly flat velocity responses, at least between 0.03 and 8.0 Hz; the sampling rates vary among 20, 40 and 100 per second. Figure 1 shows the location of the Russian bolide event, and the distribution of the 50 seismic stations, which are mainly installed in Central Asia and Southwestern Europe. Parameters of these stations are listed in Table 1. The seismograms generated by the bolide explosion contain an abundance of low-frequency content. To isolate the low-frequency components and emphasize the Rayleigh wave, we apply a bandpass filter to the data with periods between 8 and 25 s. Figure 2a illustrates the obtained three-component Rayleigh waveforms recorded at station ABKT. Because the Rayleigh wave is composed of coupled P and SV waves (Edwards et al., 2008), its maximum amplitudes appear on the vertical and radial components, and the transverse component is very weak. Shown in Figure 2b are root-mean-square (RMS) amplitudes for vertical, radial, and transverse components, all of which decrease with increasing distance.

Rayleigh-wave Magnitude Measurement
By using 20 s Rayleigh waveforms, Gutenberg (1945) built a surface-wave magnitude formula that has been widely used. To study small earthquakes, Bonner et al. (2003) measured their magnitudes by using 7 s Rayleigh waves, which largely increased the number of relevant observations at regional distances. Similarly, Taylor et al. (2003) extended measurements to short-period (between 6 and 12 s) Rayleigh waves, while Russell (2006) developed a more robust measurement, which can use waveforms with a large range of frequencies. Based on three datasets,  demonstrated that the method by Russell (2006) is capable of decreasing the scatter in the magnitude estimates. Furthermore, based on these new technologies, the surface-wave magnitudes of small events can be investigated by using shortperiod Rayleigh waves at regional distances (Bonner et al., 2008;Fan N et al., 2013).
Based on the vertical-component Rayleigh waveforms, we measure the for the 15 February 2013 Chelyabinsk bolide explosion using the time domain method given by Russell (2006). The Rayleigh-wave magnitude can be represented by where is the zero-to-peak amplitude of the Butterworth-filtered surface wave in nanometers, is the epicentral distance in degrees, is the period in seconds, and is the cutoff frequency of a two-pass third-order zero-phase Butter- worth bandpass filter with corner frequencies and . This technique, often referred as (VMAX), is employed for variable period, maximum-amplitude magnitude estimations Russell, 2006). Moreover, this method is not only capable of decreasing the scatter in the magnitude estimates, but also effective for estimating magnitudes of small events (Bonner et al., , 2008Fan N et al., 2013).
We first remove the instrument response from vertical-component waveforms and then extract Rayleigh waves with a group velocity window between 5.0 and 2.0 km/s. We then apply a fourthorder zero-phase Butterworth filter to obtain a series of Rayleigh waveforms in variable periods between 8 and 25 s. Next, we use an envelope function to measure the maximum zero-to-peak amplitudes from the filtered ground displacements, and finally the results are used in equation (1) to calculate the surface-wave magnitudes.
The processing of the vertical-component Rayleigh waveform at station WMQ is shown as an example, in Figure 3. The obtained peak amplitudes at different periods vary between 270.91 and 484.89 nm, and the magnitudes range from 3.62 to 4.41. Following Russell (2006), we choose the maximum-amplitude magnitude as the station-event magnitude. All station-event Rayleigh-

Estimating the Yield of the Bolide Explosion
A number of methods have been employed to estimate the seismological yields of underground nuclear explosions based on the surface-wave magnitude (e.g., Nuttli, 1986;Stevens and Murphy, 2001;Zhao LF et al., 2008;Bonner et al., 2008). The shock wave generated by an air burst, however, attenuates rapidly as it propagates through the atmosphere, and coupling between the atmosphere and solid ground is very weak. Therefore, only a frac-tion of the energy from the explosion is observed through the seismic waves. As the formulas proposed by the above mentioned studies were developed for underground nuclear explosions, they might not be directly suitable for the current purpose.
For the 1908 Tunguska explosion, Chyba et al. (1993) suggested that the best energy estimate was based on seismic records combining with the nuclear airbursts of comparable yield at the same height and the same location. Accordingly, we can utilize certain simulation results as a reference. Harkrider et al. (1974) used several theoretical source models to compute Rayleigh-wave magnitudes for high-altitude explosions over Earth structures in both continental and oceanic regions. They concluded that surfacewave magnitudes could be used to estimate the explosive yields   Figure 3. Calculation of the Rayleigh-wave magnitude from regional seismograms. Shown on the left are normalized Rayleigh waveforms bandpass filtered for periods centered between 8 and 25 s. In the center, the maximum amplitudes of these waveforms are measured from their envelopes. On the right are Rayleigh-wave magnitudes computed at different periods.
of atmospheric events if and only if an independent estimate of the burst height was obtained. With an independent estimate of the explosion height from NASA, we use theoretical relations of Harkrider et al. (1974) to estimate the Chelyabinsk bolide yield from the Rayleigh-wave amplitudes.
Based on theoretical simulations (Harkrider et al., 1974), Figure 4 illustrates the relationship of explosion yield versus height and Rayleigh-wave magnitude for the continental area. The thin curves are for yields and black dots are calculated values. It appears that Rayleigh wave excitations are complex for explosions below 10 km, but at altitudes above 20 km their trend becomes more regular. For the 15 February 2013 Chelyabinsk bolide explosion, given its height of 23.3 km and Rayleigh wave magnitude , the yield can be obtained (shown as the intersection of the two red lines). The small box at the intersection contains 4 nearby calculated values. We use these four numbers to calculate the yield by linear interpolation; the obtained yield is 686 kt.

Focal Mechanism
For most atmospheric meteorite explosions, the corresponding ground motions are triggered by the shock waves from the explosion rather than from the direct impact of the surviving meteorite (e.g., Edwards et al., 2008;Ceplecha and Revelle, 2005). Unlike an underground nuclear explosion, the amount of spall and the secondary tectonic effects of an atmospheric event are minimal, and its deviatoric moment tensor components are difficult to observe (Antolik et al., 2014). Therefore, the seismic records of this event are similar to waves caused by pressure changes in the fluid atmosphere over the solid Earth.
One of the most widely accepted source representations for seismic waves generated by atmospheric nuclear explosions is a point impulse (e.g., Ben-Menahem, 1975;Chyba et al., 1993;Langston, 2004). In these studies, the source modeled by a vertical point im-pulse impinging on a layered elastic Earth model was able to successfully account for various phases in the generated seismic waveforms; consequently, the pulse can be interpreted as the initial vertical momentum of the blast wave impact. For the 2013 bolide event, the absence of significant bolide signals on most transverse component records indicates that the source mechanism was horizontally isotropic and that the source area of Rayleigh waves was small compared with the area over which airborne waves could have been observed. Thus, we try to simulate the source as a concentrated single force to investigate the seismic focal mechanism.
We use the location (55.15°N, 61.41°E) provided by the Incorporated Research Institutions for Seismology (IRIS) as the epicenter of this event. Figure 5 shows the relationship among the airburst explosion, the epicenter, and the radiated seismic waves. The bolide explosion above Chelyabinsk excited seismic waves when its shock waves struck the ground. We simplify the impact from the shock waves as a single force that acted vertically downward on the ground surface at the epicenter.
We used seismic data from 28 of the 50 stations shown in Figure 1 to study the focal mechanism solution for the impact of shock waves from the bolide explosion. Synthetic seismic records were generated and sampled using the same group velocity window used for the real Rayleigh wave data. The synthetic seismograms were computed by the frequency-wavenumber (FK) synthetic seismogram package fk3.2 (Haskell, 1964;Wang CY and Herrmann, 1980;Zhu LP and Rivera, 2002). We adopted a single-force point source with a trial magnitude of 1.0×10 20 N at a depth of 0.01 km in a layered half-space specified by the Crust 1.0 model (Laske et al., 2013). Figure 6 compares the synthetic (red) and observed (black) Rayleigh waves, with their amplitudes normalized. A visual inspection reveals that the synthetic waveforms fit the observed waveforms quite well. Quantitatively, their correlation coefficients range from 0.51 to 0.91 with an average value of 0.74, which validate the single-force focal mechanism used for seismic Rayleigh waves generated by a bolide explosion.
The amplitudes of the synthetic Rayleigh waves are much larger than those from observed Rayleigh waves. To determine the actu-    6. Comparison between the observed (black) and the synthetic (red) Rayleigh wave seismograms. The station names and corresponding epicentral distances are indicated to the left of traces. The correlation coefficients between the recorded and the synthetic waveforms are listed to the right of these traces.

Discussion and Conclusions
We investigated the seismic characteristics of the 15 February 2013 Chelyabinsk bolide explosion. The Rayleigh-wave magnitude was 4.17±0.31 for this event, and the yield of this bolide explosion was estimated to be 686 kt. Since the synthetic Rayleigh waves fit the observations well, a single force was deemed an appropriate simplifying assumption to simulate the effect of shock wave impact.

Magnitude of the Chelyabinsk Bolide Explosion
The regional Rayleigh-wave magnitude method has three main advantages (Russell, 2006). First, this technique allows us to visually identify the phases by measuring the surface-wave amplitudes in the time domain. Second, it allows the surface-wave magnitudes to be measured at regional distances, where the 20 s Rayleigh waves required by the traditional magnitude measurement are often unavailable. Third, the application of narrowband Butterworth-filtering techniques appropriately handles Airy phase phenomena. Finally, the regional Rayleigh-wave magnitude technique is able to use seismic data recorded at any epicentral distance, whereas the traditional method is able to use data only at particular distances.
We also used the traditional method provided by Gutenberg (1945) to calculate the surface-wave magnitudes (listed in Table 1). These magnitudes range from 2.84 to 4.69 with an average value of 3.61±0.37. The average value computed by the regional Rayleigh-wave magnitude technique is larger than that from the traditional method by 0.68 magnitude unit, because the method utilized in this paper measures the Rayleigh-wave magnitude where the signal is the largest. The smaller standard deviation of the magnitudes derived from the regional method also verifies the high stability of this approach.
The body-wave magnitude (Lg) determined by the National Earthquake Information Center of the United States Geological Survey for this event was 4.2; in contrast, an (Lg) of was obtained by a formula applicable to underground explosions in East Kazakhstan (Nuttli, 1986). As the epicenter of the 2013 bolide event was near East Kazakhstan, the formula suggested by Nuttli (1986) is deemed more reasonable for the calculation of the body-wave magnitude. We therefore accept as the body-wave magnitude of this event. Consequently, we find that the surface-wave and body-wave magnitudes for the bolide event only marginally satisfy the empirical relation between and for Eurasia (Murphy et al., 1997). Two possible reasons for this bias can be explored. First, the characteristics of the emplacement medium and the zero source depth may contribute to a larger surface-wave magnitude (Bonner et al., 2008); second, unlike typical tectonic earthquakes or underground explosions, short-period signals above the noise level are difficult to observe. Thus, the body waves from this event were not recorded clearly even at the closest station, resulting in unusually small body-wave magnitude readings.
The surface wave amplitudes generated from this event were much larger than those generated from earthquakes or nuclear explosions with similar yields. Although the Rayleigh waves excited by the 2013 bolide event had amplitudes approximately three times larger than those produced by the North Korean nuclear underground explosion in 2013 at comparable distances in a particular frequency band (Heimann et al., 2013), we cannot affirm whether the earthquake energy of the Russian bolide event was larger than that of the North Korean nuclear test (Zhao LF et al., 2014), which may be due to the different frequency contents of the two source types.

× 10 14
Some scientists estimated the equivalent yield of the bolide event to be as low as 0.1 kiloton, while others reported that it could be nearly 57 Mt. When the area of the shock wave hit at the ground surface is assumed to be 100 km 2 , the energy of the shock wave will be equal to J (71.8 kt) (Chernogor and Rozumenko, 2013). In reality, the area was larger than 100 km 2 (e.g., Emel 'yanenko et al., 2013;Popova et al., 2013); as a result, the energy of the bolide explosion was probably greater than 71.8 kt.
When a bolide explodes, the energy carried by the ensuing shock waves can be used to estimate the explosion's total energy. Part of energy in the shock waves will transfer into seismic wave energy after the shock waves impact on the ground. Apparently, simulating the propagation of shock waves in Earth's atmosphere can play an important role in improving estimates of the explosive yield. However, we were unable to obtain detailed atmospheric data on the day of the bolide event; and even if such data were available, the attenuation of shock waves in the atmosphere is extremely difficult to determine, as the states of the atmosphere at different depths change rapidly with time. Moreover, the air to ground coupling coefficient currently available is not generic.
For these reasons, we adopted the theoretical relations among the explosion height, Rayleigh-wave magnitude, and explosion yield. Similar to the formula used herein to calculate the surfacewave magnitude, the formula in Harkrider et al. (1974) uses variable-period surface waves in the traditional format. Some disagreement between the surface-wave magnitudes calculated by these two methods exists, but this disparity is not sufficient to affect the estimated yield significantly. In the region of the intersection point located in Figure 4, the peak energies of the Rayleigh waves from this event and the simulated Rayleigh waves in Harkrider et al. (1974) are all located at periods of approximately 20 s. Upon comparing our findings with previous results, our estimation agrees with the yield (approximately 500 kt) accepted by most agencies.
However, even with these considerations and assumptions, notable uncertainties are still involved when estimating the yield. The relations among the Rayleigh-wave magnitude, explosion yield, and explosion height are modeled by both the Gutenberg continental model, which is taken from Ben- Menahem and Harkriker (1964), and an atmosphere model represented by isothermal layers composing the standard ARDC atmosphere (Wares et al., 1960). The continental model and atmospheric model cannot exactly represent the geological properties within the area of investigation or the real atmospheric conditions at the time of the bolide event, and thus some bias in the yield is introduced when directly applying these models.

Focal Mechanism
There are slight differences between the observed and synthetic Rayleigh waves. One reason for this could be that the single-force model cannot exactly represent the seismic source of the bolide event because the focal region is located across an area rather than at a single point and the shock waves do not impact on the ground at an exactly 90 degree incidence angle. Another reason for the difference may be differences between the velocity model used in the calculation and the real velocity structure. In spite of this, high correlations between the observed and synthetic Rayleigh waves reveal that a single force plays a leading role in generating waves, and thus can represent a good equivalent source for this seismic event.

Impact Pressure of the Shock Waves on Ground
Because shock waves attenuate rapidly as they propagate through the atmosphere, they impact on a limited region. The maximum distance from the epicenter to the significantly damaged regions in this event ranged between 30 km and 100 km (e.g., Emel'yanenko et al., 2013;Lobanovsky, 2014;Chernogor and Rozumenk, 2013;Popova et al., 2013). Based on a map of the extent of glass damage on the ground provided by Popova et al. (2013), we can observe that the maximum radius of damage at the surface was approximately 100 km, but the location of significant damage is within a radius of 30 km. Assuming that the earthquake was mainly caused by the shock wave pressure loaded on an area of no more than 30 km radius, and combining this assumption with the single-force model obtained above, we estimate the average pressure to have been 0.36 kPa. In contrast, the pressures on central Chelyabinsk provided by most researchers (e.g., Emel'yanenko et al., 2013;Brown et al., 2013;Chernogor and Rozumenko, 2013;Avramenko et al., 2014) range from 0.7 kPa to 3.8 kPa. Apparently, our estimation of the average pressure is smaller than the value of the previously determined pressure at central Chelyabinsk. As shown in Figure 7, one of the most important reasons for this is the rapid attenuation of shock waves with the increase of distance. Furthermore, due to the huge imped-ance contrast between the air and the solid Earth, only a small fraction of the energy will be converted from shock wave to seismic wave. Based on the pressure data, Chernogor and Rozumenko (2013) suggest that the shock waves from the bolide explosion remained strong enough to cause partial destruction within a radius of up to 100 km. This range of destruction is equal to what we have accepted; therefore, we use the pressure data provided by Chernogor and Rozumenko (2013) to estimate the impacting force of the shock waves on the ground (Figure 7). By numerically integrating the empirical pressure equation in Figure 7, we obtained a total impact force of 2.58×10 13 N . In contrast, the force that fit with the seismic data is 1.03×10 12 N, much smaller than the directly estimated shock wave pressure. We thus estimate that only 4.0% of the shock wave energy impacting on the ground was converted into seismic energy. In the empirical equation fitting the data, P and D represent the pressure and epicentral distance, respectively.