Local magnitude scale in Slovenia

In the paper a calibration study of the local magnitude scale in Slovenia is presented. The Seismology and Geology Office of the Slovenian Environment Agency routinely reports the magnitudes MLV of the earthquakes recorded by the Slovenian seismic stations. The magnitudes are computed from the maximum vertical component of the ground velocity with the magnitude equation that was derived some thirty years ago by regression analysis of the magnitudes recorded by a Wood-Anderson seismograph in Trieste and a short period seismograph in Ljubljana. In the study the present single magnitudeMLV equation is replaced by a general form of the Richter local magnitude MWA equation. The attenuation function and station-component corrections that compensate the local effects near seismic stations are determined from the synthetic Wood-Anderson seismograms of a large data set by iterative least-square method. The data set used consists of approximately 18 000 earthquakes during a period of 14 yr, each digitally recorded on up to 29 stations. The derived magnitude equation is used to make the final comparison between the newMWA magnitudes and the routinely calculated MLV magnitudes. The results show good overall accordance between both magnitude equations. The main advantage of the introduction of station-component corrections is the reduced uncertainty of the local magnitude that is assigned to a certain earthquake.


Introduction
The Seismology and Geology Office of the Slovenian Environment Agency is recording and processing data about the earthquakes occurring in Slovenia and the surrounding areas.For the last 30 yr the M LV magnitudes of the earthquakes recorded by the Slovenian seismic stations were computed using the maximum ground velocity on the vertical component (Cecić et al., 2005;Bormann at al., 2002) with the geometrical spreading and attenuation function determined to fit the Richter (1935Richter ( , 1958) ) magnitudes M WA recorded by a Wood-Anderson seismograph in Trieste.In the last 12 yr the digital seismic network of the Republic of Slovenia has grown from 6 to 28 permanent seismic stations (Fig. 1) and recorded large set of earthquakes, but attenuation relation for the local magnitude has not been re-evaluated yet.Therefore, the need for calibration or at least verification of the presently used magnitude equation is evident, as similar studies have been carried out for neighboring regions, for example in Italy (e.g.Bindi et al., 2005;Bragato and Tento, 2005).
Out of the large data set of earthquakes occurring in the area between 44.6 • and 47.0 • N latitude and 12.0 • and 17.0 • E longitude in the period from January 1997 to December 2010 a carefully selected high quality sub-set is used to calibrate the M WA magnitude equation.An iterative least-squares method is used to determine distance attenuation coefficient and station-component corrections 1 in the M WA magnitude equation from more than 33 000 automatically determined amplitudes from synthetic Wood-Anderson seismograms of more than 1800 earthquakes (Fig. 2).The magnitude span of the data used is between M WA = 0.6 and M WA = 5.2, with all but few tens of events between M WA = 1 and M WA = 3 and the hypocentral distance span from 20 to 280 km, with less than 0.4 % of the distances larger than 210 km.
Besides the data from the Slovenian stations the data from the Italian seismic station Trieste (TRI) is used in the study. 1 The term station-component correction is used to stress that each of the three components (EW, NS, or vertical) of the seismogram is usually assigned its own correction.Nevertheless a shorter term station correction will often be used as a synonym.
J. Bajc at al.: Local magnitude scale in Slovenia de has not been re-evaluated yet.Therefore, the need for calibration or at least verification resently used magnitude equation is evident, as similar studies have been carried out for ring regions, for example in Italy (e.g., Bindi et al., 2005;Bragato and Tento, 2005).In part this is because of the historical importance of the TRI station for the derivation of the local magnitude equation for Slovenia and in part in order to be able to compare the results of this study to the magnitudes, obtained in the neighboring countries.Since most of the regional seismicity is constrained to regions close to Slovenian Italian border a comparison of the magnitudes calculated by the two seismic networks seems a natural choice.
The paper proceeds with the section about the methodology and data.The criteria for data selection and the methods used to obtain the unknown parameters of the calibrated magnitude equations are described in detail.In the following section the obtained attenuation function and station corrections are presented.Finally the results are discussed and summarised in the conclusions.

Methodology and data
Currently, the M LV of an event in Slovenia or its vicinity is determined as the mean of the M LV magnitudes at individual stations, using a single magnitude equation where A is the ground displacement amplitude in nanometers, T the period in seconds and r the hypocentral distance of the station in kilometers.As the modern seismometers measure ground velocity the vertical component of the seismogram is searched for the highest peak-to-peak value that corresponds to (A/T ) max and the value of the ground displacement amplitude A is calculated from the obtained maximum ground velocity (Bormann at al., 2002).The magnitude magnitudes calculated by the two seismic networks seems a natural choice.
The paper proceeds with the section about the methodology and data.The criteria for dat tion and the methods used to obtain the unknown parameters of the calibrated magnitude eq are described in detail.In the following section the obtained attenuation function and station 45 tions are presented.Finally the results are discussed and summarised in the conclusions.

Methodology and data
Currently, the M LV of an event in Slovenia or its vicinity is determined as the mean of th magnitudes at individual stations, using a single magnitude equation where A is the ground displacement amplitude in nanometers, T the period in seconds an hypocentral distance of the station in kilometers.As the modern seismometers measure velocity the vertical component of the seismogram is searched for the highest peak-to-pea that corresponds to (A/T ) max and the value of the ground displacement amplitude A is cal from the obtained maximum ground velocity (Bormann at al., 2002).The magnitude equation equation (Eq. 1) has been developed approximately thirty years ago by Ribarič, then the head of the Seismological Survey of Slovenia, but no written documentation of the methodology and data used is available.He developed the M LV equation (Eq. 1) by adjusting the numerical constants in it in such a way to obtain for the same earthquakes the best agreement with the available M WA values from the station TRI, then equipped with a Wood-Anderson seismograph.
An appropriate ansatz for the M WA magnitude equation to replace Eq. ( 1) should originate from the Richter (1935Richter ( , 1958) ) magnitude equation as introduced for Southern California.For numerical modelling the tabulated attenuation function in the original paper (Richter, 1935) is not suitable.Therefore the analytical form of the M WA magnitude equation which is based on the original Richter idea, but was developed for Southern California several decades later by Hutton and Boore (1987) and Boore (1989), is taken as the basis of the ansatz.In Eq. (2) A is the maximum amplitude on one of the horizontal components of the Wood-Anderson seismogram in millimeters and r is the hypocentral distance in kilometers.Since there is no Wood-Anderson instrument in the data set used, all Wood-Anderson seismograms are synthetic, simulated from broad-band seismometer data.It is easily verified that the Eq. ( 2) satisfies the Richter definition of the magnitude: an earthquake that produces the amplitude A of 1 mm at the station that is 100 km away, is assigned are too many C j to show all these results for different years in this short paper, but an idea ocedure is illustrated with the determination of the parameter a. Results for different subummarised in Table 1.There is no a priori reason that a, which is related to geometrical and attenuation of the seismic waves in a region, would change in time.So, the change alue, obtained from data sets within different time intervals, is related to the uncertainty of od and of the input data.Therefore, the fluctuations of the a values in Table 1 provide an e uncertainty of the value of the parameter a ± 0.04 = 1.38(1 ± 0.03) .
(7) 7 a magnitude M = 3.0.The constants 1.110 and 0.00189 are obtained by fitting the attenuation function to Southern California data (Hutton and Boore, 1987;Boore, 1989).Of course the attenuation in Slovenia is different from the one in Southern California, but the ansatz for the M WA magnitude equation for j -th station-component with added station correction C j should still be in agreement with Eq. ( 2), therefore it reads where constant 3.0 is retained for easier comparison of the original Richter (Eq.2) and the proposed (Eq. 3) magnitude equation.Note that since there are two horizontal components on each station there are two times more stationcomponent corrections than stations.In Eq. (3) a and b depend on the regional attenuation and geometrical spreading of seismic waves.
In order to keep the explanation of the methodology general, let us denote the individual magnitude of the i-th earthquake on the j -th station as M ij and omit the indices WA for now.We define the magnitude of the i-th earthquake M i as the mean value of n available individual magnitudes where n varies from earthquake to earthquake depending on its location, magnitude, and number of stations operating at the time.We name this magnitude the earthquake magnitude in order to distinguish it from the individual magnitudes of the same earthquake.
Station corrections C j are introduced to reduce the systematic differences between individual magnitudes of the same earthquake due to local effects related to each seismic station.That is why the station corrections are determined from the condition that the sum of squared differences between individual magnitude and earthquake magnitude for  6).In the first column the time period is indicated, in the second one the number of amplitude data points in the time period is indicated, in the third column the number of the earthquakes providing amplitude data is written, and finally in the fourth column the corresponding value of a is written.Italic numbers in the third row are calculated from the entire data set and represent the end results of the study.

Time
No   all earthquakes is minimal where M i is calculated from Eq. ( 4).The unknowns of the problem posed are attenuation coefficients a and b, all station corrections C j , and all earthquake magnitudes M i .For a given data set of approximately 30 stations and several thousands of earthquakes the majority of the unknowns are the earthquake magnitudes.This makes the problem numerically difficult to manage and we used a simplification to overcome the large number of unknowns.Before tackling the problem further, one important point has to be addressed.Using the ansatz (Eq. 3) when minimising the sum of squared differences (Eq.5) between individual (M ij ) and earthquake (M i ) magnitudes results in a www.adv-geosci.net/34/23/2013/Adv.Geosci., 34, 23-28, 2013 solution that is indeterminable up to an additive constant C. As this constant is added to all individual magnitudes, its value is transferred from individual magnitudes to the earthquake magnitude unchanged, leaving the difference between any M ij and M i unaffected.It fallows that there is an infinite number of sets of C j values that equally well fulfil the minimal condition (Eq.5): a set of C j values produces exactly the same sum (Eq.5) as the set C j = C j + C, i.e. the set C j and the set C j are both either a solution of the problem or are both not the solution of the problem.If no additional constraint is imposed, solution may become numerically unstable -a numerical drift of the C j values during the iterative solving of the problem may occure.To overcome this numerical instability some sort of constraint has to be applied to station corrections C j .Two standard approaches are at hand.One is to fix C j for a particular station, so all the other C j 's are calculated relative to the fixed one.The other plausible possibility is to fix the average of all the station-component corrections to a fixed value, most usual choice being zero.
We follow the latter of the two approaches.
The problem is solved in the least-squares sense, numerically employing an iterative singular value decomposition method when manipulating system of equations.As the iterative method is used, there is a relatively simple way of reducing the number of unknowns.During each iteration the sum in Eq. ( 5) should be minimised by adequately changing the unknowns: a, b, C j , and M i .If the values in a particular iteration step are close enough to the final numerical solution, the earthquake magnitude values M i from previous iteration are not very different from the "real" magnitudes that are to be determined.Under this assumption the M i values in k-th iteration may be taken as constants in Eq. ( 5), having the values, obtained in (k − 1)th iteration.This reduces the number of the unknowns enormously, from several thousands to only around 60, because around 30 stations are in the data set used.Several tests were done on smaller samples to confirm that the full inversion and the proposed simplified one give exactly the same results, but the latter is faster.The minimisation is therefore done in two steps for each iteration.First the earthquake magnitudes from previous iteration (M i (k − 1)) are used as constants in Eq. ( 5) while new station correction values C j (k) are determined by minimising the sum in Eq. ( 5).To prevent the numerical drift, mentioned earlier, all C j (k) are shifted for a constant δ(k), determined from the condition that the average of all values C j (k) (C j (k) = C j (k)+δ(k)) remains zero in each iteration.In the following step new values of station corrections C j (k) are used to re-calculate earthquake magnitudes M i (k) using Eq. ( 4).The procedure is repeated until the values of the parameters a, b, and C j do not change significantly any more.As a by-product the earthquake magnitudes M i are obtained as well, but this is not that important, because they can be calculated for any event independently with obtained a, b, and C j using Eqs.( 3) and (4).The main improvement in the proposed M WA earthquake magnitude determination is through the station-component corrections in the magnitude equation (Eq.6).The uncerta magnitudes, assigned to a particular earthquake, is reduced on average for approximately Due to large initial data set of several 100 000 seismograms the only way to manage the data is automatic reading of maximum amplitudes.Automatic reading from synthetic horizontal Wood-Anderson seismograms resulted in nearly 65 000 values of A for M WA calculation.For the sake of quality control the period T of each swing with the maximal amplitude is stored as well.Because of the automated reading the quality control of the obtained A and T data is crucial.After several tests the adopted criteria of minimal number of stations, minimal signal to noise ratio, particular period interval of maximal phases, and minimal distance between the earthquake and the station resulted in 33 165 amplitude data from 1852 earthquakes.
In short, the selection process is done in the following way.First the events that are to close in time -the time window used is 30 s -are removed from the data set, because the seismograms of two such events could overlap at some stations.Next the events with calculated depths over 25 km are removed, because such events are very rare in broader Slovenia region, therefore large depths may quite often be artefacts, in particular, if an earthquake occurs outside the network.Next the minimal hypocentral distance of 20 km is selected and all the data from stations closer to an earthquake are removed.There are two reasons for the introduction of the minimal hypocentral distance.Firstly, the depths of the events carry the largest uncertainty and for small hypocentral distances the impact of depth uncertainty on the hypocentral distance uncertainty is large.Secondly, the hypocentral distance r enters the magnitude ansatz (Eq. 3) in log function so the same absolute value of uncertainty results in a much higher uncertainty of the log(r) for smaller r when compared to larger r.Next the amplitude data are filtered according to the period T .Since only the local events (up to at most a few hundreds kilometers) are considered, data with periods larger than 1.0 s Adv.Geosci., 34, 23-28, 2013 www.adv-geosci.net/34/23/2013/are removed, as longer periods may indicate artefacts.At the other end data with the periods smaller than 0.03 s are also removed, because the numerical rounding error when writing periods in the input data with two digits and having sampling rate of the order of a hundred samples per second may produce up to 50 % error.Finally, in order to reduce artefacts still further and increase the quality of the retained data, signal to noise ratio of the data used in the inversion is set to at least 6.0 and the minimal number of amplitude data per event is set to 6, regardless of the number of participating stations.
Since the distance interval of the data is mostly from 20 to 200 km (Fig. 3) with the average of approximately 76 km, the use of the two parameters a and b in the M WA ansatz (Eq. 3) is not significantly contributing to the reduction of the final root-mean-square (RMS) error.The two parameters appear to be more or less anti-correlated -increase of one results in decrease of the other, but the fit remains almost unchanged for the distance interval between 20 and 200 km.The M(r) function at smaller distances (r ∼ 100 km) is changing primarily because of the log(r/100) term, while the linear term (r − 100) prevails at large distances (r 100 km).Since the parameter a is related to M ∝ log(r/100) and b to M ∝ (r − 100) the parameter b is set to b = 0 and the general ansatz (Eq. 3) is simplified to read The ansatz (Eq.6) is the one actually used in the inversion and discussed in the rest of the paper.

Results
In order to check the uncertainty as well as the stability of the solution, several sub-sets of data are used.Typically the inversion is performed on the data for a particular year and the results for different years are compared.There are too many C j to show all these results for different years in this short paper, but an idea of the procedure is illustrated with the determination of the parameter a. Results for different sub-sets are summarised in Table 1.There is no a priori reason that a, which is related to geometrical spreading and attenuation of the seismic waves in a region, would change in time.So, the change in the a value, obtained from data sets within different time intervals, is related to the uncertainty of the method and of the input data.Therefore, the fluctuations of the a values in Table 1 provide an idea of the uncertainty of the value of the parameter a a = 1.38 ± 0.04 = 1.38(1 ± 0.03) . (7) The results of the inversion for the station-component corrections C j obtained from the data of the entire period 1997-2010 are summarised in Fig. 4. Similarly to the determination of the uncertainty of the a value the uncertainty of the station corrections are estimated.In an ideal case the stationcomponent corrections C j would completely compensate the local effects of each station-component amplitude readings and each individual magnitude M ij would be equal to the earthquake magnitude M i .Once the station corrections C j are determined, the size of the differences between individual magnitudes and earthquake magnitudes provides a measure of how reliable an individual magnitude calculated from a single amplitude reading is.As a measure of the uncertainty the uncertainty σ M of the individual magnitude M ij at the j -th station-component is defined as where n is the number of all the amplitude data at the j -th station-component, used in the inversion.In In order to compare both magnitude equations the earthquakes for which there was enough amplitude data available on vertical ground velocity seismograms as well as on simulated horizontal displacement seismograms are selected.Next both magnitudes are calculated for all 922 such events in the time interval of 14 yr, from 1997 to 2010.For each earthquake both magnitudes are drawn as a point coordinates in Fig. 5. Linear regression is used to find the relation between the two magnitudes.Relatively high value of the correlation coefficient (R 2 = 0.941) indicates that both magnitudes are strongly correlated M WA = 1.064MLV − 0.172 (10) M LV = 0.885M WA + 0.257 .
As the slope of the linear function relating M WA to M LV is nearly 1, it seems that the thirty years old original M LV can be transcribed to M WA almost one-to-one.In order to check this further, the ansatz is tested.It turns out that M = −0.057minimises the RMS error, indicating that the two magnitudes are indeed nearly the same.Even more, the ansatz (Eq.12) with the minimizing value of M = −0.057has the RMS error only 2.9 % larger than the linear regresion solution (Eq.10).The values of both magnitudes are on average almost the same, but there is an important difference: the uncertainty σ WA of the calculated M WA is significantly smaller then the uncertainty σ LV of the calculated M LV .
Contemporary single C M LV magnitude equation and new multiple C j Wood-Anderson M WA magnitude equation with suitably selected offset value M give almost equal earthquake magnitudes.This means that the attenuation function for the seismic waves, originating from local earthquakes, in Slovenia is well described by the 1.52 log(r) dependence in the presently used M LV magnitude equation (Eq.1).This is rather surprising, because of the fact that this dependence has been derived by the comparison of the data for only two seismic stations (TRI and LJU) and only small number of earthquakes.Although the attenuation function in the proposed M WA magnitude equation is described by the 1.38 log(r) dependence, which seems significantly different, the obtained M WA magnitudes do not differ significantly from the M LV values.This may in part be a consequence of the fact that M LV is based on the vertical velocity amplitudes, whereas the M WA is based on the horizontal displacement amplitudes.
The main improvement in the proposed M WA earthquake magnitude determination is achieved through the stationcomponent corrections in the magnitude equation (Eq.6).The uncertainty of the magnitudes, assigned to a particular earthquake, is reduced on average for approximately 35 %.

Fig. 1 .
Fig. 1.The stations (open squares) of the Seismic network of the Republic of Slovenia that are used in the study.Black squares indicate seismic stations Trieste (TRI) and Ljubljana (LJU) that were used in the determination of the presently used M LV magnitude equation (Eq.1).

Fig. 2 .
Fig. 2. Black dots denote earthquakes from the period 1997 to 2010 that are considered in the stu earthquakes that are actually used in the calculation of the parameters in the MWA magnitude equa coloured red.

55Fig. 2 .
Fig. 2. Black dots denote earthquakes from the period 1997 to 2010 that are considered in the study.The earthquakes that are actually used in the calculation of the parameters in the M WA magnitude equation are coloured red.
istribution of the hypocentral distances r for the data, used in the inversion.The average distance is by the arrow.thedistance interval of the data is mostly from 20 km to 200 km (Fig.3) with the average of ately 76 km, the use of the two parameters a and b in the M WA ansatz (Eq. 3) is not signifntributing to the reduction of the final root-mean-square (RMS) error.The two parameters be more or less anti-correlated -increase of one results in decrease of the other, but the fit lmost unchanged for the distance interval between 20 km and 200 km.The M (r) function r distances (r ∼ 100 km) is changing primarily because of the log(r/100) term, while the m (r − 100) prevails at large distances (r 100 km).Since the parameter a is related to (r/100) and b to M ∝ (r −100) the parameter b is set to b = 0 and the general ansatz (Eq. 3) ed to read log(A) + alog r 100 km + 3.0 + C j .(6)tz(Eq.6) is the one actually used in the inversion and discussed in the rest of the paper.ltstocheck the uncertainty as well as the stability of the solution, several sub-sets of dataTypically the inversion is performed on the data for a particular year and the results for years are compared.

Fig. 3 .
Fig. 3. Distribution of the hypocentral distances r for the data, used in the inversion.The average distance is indicated by the arrow.

Fig. 4 .
Fig. 4. The station-component corrections C j for all stations (red: EW; green: NS).For the boreho ments (GOLS, LEGS, PDKS), which do not have standard EW/NS orientations, corrections for comp and 2 are shown as EW and NS respectively.Black intervals indicate the ±σ M j uncertainty intervals.

Fig. 4 .
Fig. 4. The station-component corrections C j for all stations (red: EW; NS).For borehole instruments (GOLS, LEGS, PDKS), which do not have standard EW/NS orientations, corrections for components 1 and 2 are shown as EW and NS respectively.Black intervals indicate the ±σ Mj uncertainty intervals.

Fig. 5 .
Fig. 5. Comparison of the new M WA and original M LV magnitudes for the earthquakes, used in the

235 10Fig. 5 .
Fig. 5. Comparison of the new M WA and original M LV magnitudes for the earthquakes, used in the inversion.
Fig. 4 black intervals indicate the uncertainty σ Mj for each stationcomponent correction C j .As the last step an overall comparison of the new M WA and original M LV magnitude equation is done.The new M WA magnitude equation takes the form M j WA = log (A) + 1.38 log r 100 km + 3.0 + C j .

Table 1 .
Different sub-set results for parameter a in Eq. (