Determining the focal mechanisms of the events in the Carpathian region of Ukraine

The modification of the matrix method for constructing the displacement field on the free surface of an anisotropic layered medium is presented. The source of seismic waves is modelled by a randomly oriented force and seismic tensor. A trial and error method is presented for solving the inverse problem of determining parameters of the earthquake source. A number of analytical and numerical approaches to determining the earthquake source parameters, based on the direct problem solutions, are proposed. The focal mechanisms for the events in the Carpathian region of Ukraine are determined by the graphical method. The theory of determination of the angles of orientation of the fault plane and the earthquake’s focal mechanism are presented. The focal mechanisms obtained by two different methods are compared.


INTRODUCTION
The main data sources in seismology are the seismic records of natural or man-made events that are received on the Earth surface.The task of modern seismic analysis is to obtain the maximum possible information about the nature of wave-fields propagation.Solving these problems involves the study of seismic regions of Ukraine and interpretation of wave fields in order to determine the earthquake focal mechanisms.
In recent years one of the most important methods is the development of approaches for constructing the theoretical seismograms, which allow the study of the structure of the medium and determination of the earthquake source parameters.The effects on the wave field and seismic waves' propagation in the Earth's interior should be considered when calculating these seismograms.Thus, the displacement field, which is registered on the free surface of an inhomogeneous medium, depends on the model of the geological structure and the physical processes in the source.
Interpretation of seismic research can predict the dynamic properties of elastic media, and consider the effects of anisotropy in the inversion problems of determining the source parameters.Therefore, the problem of mathematical modelling of seismic wave propagation in anisotropic medium is relevant.Over the past decade the considerable experience in theoretical and algorithmic solutions of a wide range of dynamic seismology problems is accumulated.There are plenty of methods for solving such problems, which are quite effectively used in geophysics, including seismology.Analytical problem solving methods are developed only for a relatively narrow range of tasks.More precise, and hence more complex mathematical models are implemented by numerical methods.The last give a solution only in certain limited areas of model medium, and this is the main drawback of numerical methods.This means that the use of numerical methods, including finite difference method [Fuchs 1977, Ilan 1975, Bullen 1953, Yang 2002, Zahradnik 1975] and finite element method [Thomson 1950, Woodhouse 1978] for modelling of seismic wave propagation in inhomogeneous anisotropic media gives very high accuracy results, but requires a grid which covers the entire area occupied by the investigated object and a significant amount of computer resources for the solution of highdimensional systems of algebraic equations.Therefore, it is difficult to implement, even with the use of modern computational tools, including clusters.The matrix method is used to obtain solutions, which avoid the complicated procedures to satisfy all boundary conditions.The usefulness of solutions obtained by this method is considered in [Babuska 1981, Bachman 1979, Backus 1962, Behrens 1967, Dunkin 1965].The matrix method allows for a common approach to examine the propagation of waves in a wide class of systems.This method allows to obtain solutions in a more compact and convenient form for further analytical and numerical calculations.
In the 50's of 20th century Thomson and Haskell first proposed a method for constructing interference fields by simulation of elastic waves in layered isotropic half-space with planar boundaries [Haskell 1953].The matrix method was developed in the works [Malytskyy 1998, 2008, 2010, Kennett 1972, Cerveny 2001,Chapman 2004].The stable algorithms of seismograms calculation for all angles of seismic wave's propagation are obtained.The matrix method is generalized for low-frequency waves in inhomogeneous elastic concentric cylindrical and spherical layers surrounded by an elastic medium.The concept of the characteristic matrix determined by physical parameters of the environment is developed.The matrix method is used for seismic waves propagation in elastic, liquid and thermoelastic media.In addition, it has been generalized for the study of other processes described by linear equations.The advantage of the matrix method is the ability to compactly write matrix expressions that are useful both in analytical studies and numerical calculations.
The matrix method and its modifications are used to simulate the seismic waves` propagation in isotropic and anisotropic media.This method is quite comfortable and has several advantages over other approaches.Both advantages and disadvantages of the matrix method are well described in [Malytskyy 2010, Thomson 1950, Ursin 1983, Thomsen 1966].
Today in seismology much attention is given to mathematical modelling as one of the main tools for the analysis and interpretation of the wave fields.In this paper using a modification of the matrix method Thomson-Haskell, rigorous equations for the wave field on the free surface of inhomogeneous anisotropic medium are obtained, when a source of seismic waves is located within a homogeneous anisotropic layer and presented by the seismic moment tensor.Note that the problem of wave fields modelling generated by a source, which is presented in terms of seismic moment tensor, also has practical applications in seismology.Using this method, the approaches to determining the displacement field are developed for different types of earthquake sources that will be shown in the following sections.

DIRECT PROBLEM
The problem of wave fields modelling, when the source is presented by seismic tensor moment, has practical applications in seismology.Therefore, the development of methods for determining the displacement field on the free surface of an anisotropic inhomogeneous medium for sources of this type is an actual task and needs to be resolved.
In this section the propagation of seismic waves in inhomogeneous anisotropic medium is considered.The modification of the matrix method of construction of wavefield on the free surface of an anisotropic medium is presented.The earthquake source represented by a randomly oriented force or a seismic moment tensor is placed on an arbitrary boundary of a layered anisotropic medium.The theory of the matrix propagator in a homogeneous anisotropic medium by introducing a "wave propagator" is presented.It is shown that for anisotropic layered medium the matrix propagator can be represented by a "wave propagator" in each layer.The displacement field on the free surface of an anisotropic medium is obtained from the received system of equations considering the radiation condition and that the free surface is stressless.

Theory of modification of the matrix method
The problem of wave fields modelling, when the source is presented by seismic moment, has practical applications in seismology.Therefore, the development of methods for determining the displacement field on the free surface of an anisotropic inhomogeneous medium for sources of this type is an actual task and needs to be resolved.
In this paper the propagation of seismic waves in anisotropic inhomogeneous medium is modelled by system of homogeneous anisotropic layers, as shown in (Fig. 1).The each layer is characterized by the propagation velocity of P-and S-wave and density.At the boundaries between layers hard contact condition is met, except for the border, where the source of seismic waves is located.

Fig. 1. Model vertically inhomogeneous medium
The earthquake source is modelled by nine pairs of forces, which represented a seismic moment tensor.This description of the point source is sufficiently known and effective for simulation of seismic waves in layered half-space [Haskell N.A., 1953].].In general, the source is also assumed to be distributed over time, i.e. seismic moment M 0 (t) is a function of time.This means that the physical process in the source does not occur instantaneously, but within a certain time frame.It is known for our seismic events (Mw ~ 2-3) that the time during which occurred the event may be 0.1 -0.7 seconds.The determination of the source time function is an important seismic problem.In this chapter the direct problem solution is shown, when a point source is located on an arbitrary boundary of layered anisotropic media.
We assume the usual linear relationship between stress τ ij and strain e kl where u=(u x , u y , u z ) T is displacement vector.
The equation of motion for an elastic homogeneous anisotropic medium, in the absence of body forces is [Fryer et al, 1984] where ρ is the uniform mass density, and ijkl c are the elements of the uniform elastic coefficient tensor.
Taking the Fourier transform of (1) and (2), we obtain the matrix equation [Fryer et al, 1987] ) is the vector of displacements and scaled tractions, . With the definition of b the system matrix A has the structure ; where T, S and C are 3×3 sub matrices, C and S are symmetric.
For any vertically stratified medium, the differential system (3) can be solved subject to specified boundary conditions to obtain the response vector b at any desired depth.If the response at depth z 0 is b(z 0 ), the response at depth z is (4) where P(z, z 0 ) is the stress-displacement propagator.
To find this propagator, it is necessary to find the eigenvalues (vertical slownesses), the eigenvector matrix D, and its inverse D -1 [Fryer et al, 1984]: (5) where Q is the "wave" propagator [Fryer et al, 1984]: In the isotropic case the eigenvector matrix D known analytically, so the construction of the propagator is straightforward.In the anisotropic case, analytic solutions have been found only for simple symmetries so in general, solutions will be found numerically.
The layered anisotropic medium, which consists of n homogeneous anisotropic layers on an anisotropic halfspace (n +1) (Fig. 1), is considered.The source in the form of a jump in the displacementstress is placed on the s-boundary (Fig. 1); it is easy to write the following matrix equation, using (5-6): , where where . Using ( 7) and the radiation condition (with a halfspace (n+1) the waves are not returned), and also the fact that the tension on the free surface equals to zero, we obtain a system of equations: Using only the homogeneous equations is sufficient to get the displacement field on a free surface: The stress-displacement discontinuity is determined via the seismic in matrix form [Fryer G.J. et al, 1984]: where M xx , M yy , M zz , M xz , M yz , M yx , M xy , M zy , M zx -components of the seismic moment tensor, and c 13 , c 23 , c 33 , c 44 , c 55 -components of the stiffness matrix.
As a result, the displacement field of the free surface of an anisotropic medium is in the spectral domain as: where Using (9) and three-dimensional Fourier transform, we obtain a direct problem solution for the displacement field of the free surface of an anisotropic medium in the time domain as: where z R -epicentral distance, p x , p y -horizontal slowness.

INVERSE PROBLEM
It is known that inverse problems are inherently incorrect.In seismology methods and approaches often are used, which are reduced to the selection of the physical characteristics of the studied environment or earthquake [Вrace 1978, Clinton 2006, Cohn 1982, Hartzell 1979, Santosa 1986, Hanks 1979, Honda 1957, 1962, Scholz 1973].The development of new methods and algorithms for determining of source parameters is relevant and important issue.Of course, there is not general and reliable approach.
Furthermore, it is impossible to consider all effects in modelling wave processes during propagation of seismic waves in heterogeneous environments.Therefore, for an anisotropic medium it is difficult to construct a theory that would be based only on analytical expressions.Thus, we must resort to numerical solution of equations.

Traditional graphic method and trial and error method for determining of the earthquake source parameters
The focal mechanism solution for earthquakes in the region of low seismic activity today is the actual problem.Particularly it is very important for Carpathian region of Ukraine, where insufficient number of stations is in addition to low seismic activity.It is impossible to determine a focal mechanism by software packages.
The software packages Matlab 7.12.0(R2011a) is used for programming in this paper.Algorithms and software of the direct dynamic problem solving is based on the method described in Section 2.1.This method is based on the use of matrix and wave propagators, which are applied to inhomogeneous anisotropic media modelled by bundle of homogeneous anisotropic layers with parallel boundaries.Algorithm for calculating the displacement field on free surface of a layered anisotropic medium in the Cartesian coordinate system ( ) is based on certain physical and mathematical constraints: 1) heterogeneous anisotropic medium is modelled by a bundle of homogeneous anisotropic layers; 2) homogeneous anisotropic layers are separated by parallel boundaries; 3) contact between the layers is considered hard (continuity of displacements and stresses); 4) a point source is located within any homogeneous anisotropic layer.Thus, the algorithm for calculating the displacement field on free surface of a layered anisotropic medium defined by (10) in the Cartesian coordinate system.The equations for the direct dynamic problem are obtained by means of numerical calculations.
In the algorithms and programs is used Fast Fourier Transform to the variables (t,ω).Maximum frequency, step frequency and sampling time were chosen from the conditions: , where ; 2 Waves are excited by a point source, represented by seismic moment tensor.The relationship between the components of the seismic moment tensor and fault plane orientation angles is given as [Aki 2002]:  where М 0 =µAu(τ) -seismic moment; δ -a dip angle; φ s -a strike angle; λ -a slip angle.
The tensor ( 11) is defined by the geometric orientation of the fault plane and the value of seismic moment М 0 .
Obtaining of the analytical expressions to determine the earthquake source parameters, when the source is represented by seismic moment tensor, is difficult.The most accurate results of the inverse problem solving for the source parameters are obtained by the trial and error method.In this method, the synthetic seismograms are calculated many times for all possible combinations of orientation angles of the fault plane and the velocity model for the Carpathian region of Ukraine.The correlation coefficients are calculated for the all these synthetic seismograms and real record of event.The biggest correlation coefficient corresponds to the most probable combination of orientation angles of the fault plane.The best results of solving are for the records from stations in smaller epicenter distance, these records usually have the lowest noise level.The results obtained by this method are compared with results for the same event but received by graphical method [Malytskyy 2013a].
In the trial and error method the matrix equation ( 9) is solved for the velocity model for the Carpathian region of Ukraine (Table 1) and for the stress-displacement discontinuity (8), where components of seismic tensor are determined via oriental angles of the fault plane (11).
Table 1 The traditional graphical method based on the first arrival P-waves [Malytskyy 2013a, Baumbach 2011] using information about fuzzy first motion [Cronin 2004] and the S/P amplitude ratio [Hardebeck 2003].
The polarities first motion P-waves was defined from complete records seismograms taking into account the possible inversion of the sign on the z-component.A logarithm of the amplitude ratio S/P is calculated using data from the three components seismic records of this event at each station [Hardebeck 2003, de Natale 1994].Input data for the azimuth and take-off angle are calculated by software packages for each event.
Most often an approach is used where nodal planes are plotted on a lower-hemisphere stereographic projection such as to best fit the polarities of first arrivals of P waves at the stations location of a station polarity on the projection depending on the station azimuth and take-off angle of the ray of first arrival connecting the source and the station.
These focal mechanisms are determined using a method that attempts to find the best fit to the direction of P-wave first motions observed at each station.For a double-couple source mechanism (or only shear motion on the fault plane), the compression first-motions should lie only in the quadrant containing the tension axis, and the dilatation first-motions should lie only in the quadrant containing the pressure axis.Accuracy focal mechanism solution depends on the input data: velocity model and coordinate of the hypocenter (they determine the take-off angle), quality of seismic records and sign inversion on the seismometer, so that "up" is "down" (they determine character entry wave).S/P amplitude ratios are applicable because of P-wave amplitude being the largest on P and T axes of focal mechanism and the smallest near the nodal planes, while the S-wave amplitude being the largest near the nodal planes.S/P amplitude ratios with a wide range of values can more accurately constrain the location of seismic station projections on the focal sphere.The larger S/P amplitude ratios, the closer the location of the seismic station projection to the nodal line.
The seismic moment is computed according to: where r -hypocentral distance, p v -P-wave velocity, ρ -density, 0 u -low-frequency level (plateau) of the displacement spectrum, Θ -average radiation pattern and a S surface amplification for P waves.The source radius R is computed from the relationship: where c f -corner frequency of the P wave.
The size of the circular rupture plane is computed as: The average source dislocation is according to ) where the shear modulus computed by The stress drop, seismic energy and magnitude ML are computed according to:

Determining the parameters of earthquake sources
In approving the proposed trial and error method for determining of the earthquake source parameters the four events in the Carpathian region of Ukraine were considered.For each of these events an earthquake focal mechanism is determined by the trial and error method and graphic method.These focal mechanisms are compared.
Location map of the projection of seismic stations in the Carpathian region of Ukraine and specified epicenter of the four events in the Carpathian region of Ukraine is shown on Fig 2. Spectral parameters for the considered events in the Carpathian region of Ukraine calculated by (12)(13)(14)(15)(16)(17)(18)(19) are shown in Table 2.
The first step of solving the inverse problem of seismology is determining the parameters of earthquake source by graphic method.The data for the determining the focal mechanism by the traditional graphic method are sign of first arrival of P-wave, azimuth from epicenter on the seismic station, take-off angle and the S/P amplitude ratio on each seismic record of this event.The input data for graphic method are given in Table 3-6 for all considered seismic events.Taking into account the input data and S/P amplitude ratio (Table 3) the most probable focal mechanism of event 2012.01.06 is chosen and shown on   The most probable focal mechanism of event 2012.10.24 is determined by the graphic method using input data (Table 5) and shown on fig 5. Fig. 5. a -location of the projection of seismic stations and the nodal planes according to the input data (Table 5), b -focal mechanism of event 2012.10.24 determined by the graphic method Taking into account input data (Table 6) the most probable focal mechanism is chosen by graphic method and shown on fig.6. Fig. 6. a -location of the projection of seismic stations and nodal planes according to the input data (Table 6), b -focal mechanism of event 2013.04.04 determined by the graphic method The parameters of the focal mechanisms of the events in the Carpathian region of Ukraine determined by the graphic method are done in table 7. The second step of solving the inverse problem of seismology is determining the parameters of earthquake source by the trial and error method.Using the velocity model for the Carpathian region of Ukraine (Table 1), the wave-field on a free surface ( 10) is calculated many times for all combination of dip (δ), strike (φ s ) and slip (λ) angles.All of these synthetic seismograms are compared with real records of this event.The Pearson's correlation coefficients are calculated for all of these synthetic waveforms and real seismograms on each seismic station which record this event.The largest coefficient corresponds to the most probable combination of fault plane orientation angles (φ s , δ, λ).On Fig. 7   The focal mechanism of event 2012.01.10 is also determined by described above the trial and error method.Similarly to the previous case the synthetic seismograms and correlation coefficients are calculated.On Fig. 9 the coefficients R between synthetic and real seismograms are shown.On the horizontal axis a nodal plane identifier (combinations of dip, strike and slip angles) is plotted.The maximum correlation coefficient R corresponds to the focal mechanism.On Fig. 10. the focal mechanisms determined by the trial and error method using real seismograms from different seismic stations and average variant of the focal mechanism of event 2012.01.10 are shown.The third step is comparative analysis between the focal mechanisms obtained by different methods.Comparing the focal mechanisms determined by graphic method and trial and error method of the event 2012.01.06 (Fig 3.b,8.e);the event 2012.01.10 (Fig. 4.b, 10.e); the event 2012.01.24 (Fig. 5.b, 12.e) and the event 2013.04.04 (Fig. 6.b, 14.e), we can conclude that the results obtained by different methods are similar.Thereby, these two methods can be used for determining the focal mechanisms of the local events.

CONCLUSION
The results of this paper contribute to the fundamental understanding of wave propagation in anisotropic media.A numerical technique for computing synthetic seismograms has been developed in the framework of the theory.Wave propagation in multilayered media requires that displacement and stress vectors be continuous everywhere, including the interfaces.
Seismologists have been able to invert the rupture process of a number of earthquakes and many of the features predicted by simple dynamic source models have been quantified and observed.Foremost among these is the shape of the FF spectrum, the basic scaling laws relating particle velocity and acceleration to properties of the fault, such as size, stress drop and rupture velocity.Recent inversions of earthquake slip distributions using kinematic source models have found very complex source distributions that require an extensive reappraisal of classical source models.
It is shown that the developed method for determining the earthquake parameters can be used successfully using real records.It should be also noted that the proposed method for determining the seismic moment tensor can be used in seismology for a class of problems, when the velocity model of the medium is known.Thus, the methods, approaches, algorithms, software for the propagation of seismic waves and results of direct and inverse dynamic problems of seismology proposed and developed by the authors and highlighted in the paper, can be successfully used in the study of the seismic regions and effective implementation in the construction of the earthquake source mechanism which is crucial for seismic regions of the country.

Fig. 2 .
Fig. 2. Location map of the projection of seismic stations in the Carpathian region of Ukraine and specified epicenter of the four events in the Carpathian region of Ukraine Table 2. Spectral parameters for the events in the Carpathian region of Ukraine calculated by (12-19) event M 0 , N•m f cp , Hz R, m A, m 2 , 10 5 Fig 3.

Fig. 3
Fig. 3. a -location of the projection of seismic stations and nodal planes according to the input data (Table 3), b -focal mechanism of event 2012.01.06 determined by the graphic method

Fig. 4
Fig. 4. a -location of the projection of seismic stations and nodal planes according to the input data (Table 4), b -focal mechanism of event 2012.01.10 determined by the graphic method.
the Pearson's correlation coefficients between synthetic waveforms and real seismograms on seismic stations (BERU, BRIU, KORU, MUKU, MEZ, NSLU, and RAKU) are shown.On the horizontal axis a nodal plane identifier (combinations of dip, strike and slip angles) is plotted.The maximum coefficient R corresponds to the most probable focal mechanism.On Fig. 8. the focal mechanisms determined by the trial and error method using real seismograms from different seismic stations and average variant of the focal mechanism of event 2012.01.06 are shown.

Fig. 7 .
Fig. 7.The correlation coefficients for the event which took place near NNP "Synevyr" on 2012.01.06

Fig. 9 .
Fig. 9.The Pearson's correlation coefficients between synthetic and real seismograms for the event which took place near NNP "Synevyr" on 2012.01.10

Fig. 11 .
Fig. 11.The Pearson's correlation coefficients between synthetic and real seismograms for the event which took place near village Ugla on 2012.10.24

Fig. 13 .
Fig. 13.The Pearson's correlation coefficients between synthetic and real seismograms for the event which took place near village Nyzhnje Selyshche on 2013.04.04 . Velocity model for the Carpathian region of Ukraine

Table 6 .
Input data for the determining the focal mechanism of event 2013.04.04

Table 7 .
Parameters of the focal mechanisms determined by the graphic method

Table 8 .
Parameters of the average variant s of the focal mechanisms determined by the trial and error method