Calibration of the local magnitude scale (Ml) for Central Southern Africa

The Central Southern Africa is an area of moderate seismic activity generally caused by the presence of the East Africa Rift System. To improve the quantification of seismicity in this region, we propose a local magnitude scale (Ml), based on the original Richter definition to be used by the Zimbabwe and Botswana national networks. The magnitude scale is developed using 621 seismic events that occurred between 1997–2000 and 2013–2018. These events were recorded by 61 broadband seismic stations located in Botswana, Zimbabwe, and South Africa. We evaluated 5306 traces of zero to peak maximum amplitude, recorded on the horizontal channel from simulated Wood-Anderson seismograms. All 5306 maximum amplitude measurements were inverted simultaneously to determine the attenuation constants. The resultant Ml is defined as Ml=log10A+0.80×log10(R)+0.00086×R−1.37±S\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$ M_{l} = \log _{10}A+0.80 \times \log _{10}(R)+0.00086 \times R-1.37 \pm S $\end{document} in which A is ground displacement amplitude determined from instrument corrected synthetic Wood-Anderson seismograms in nanometres, R is the epicentral distance (km), and S is the station correction. Station corrections were determined for all stations during analysis resulting in values ranging between − 0.44 and + 0.31. The range in station corrections suggests a strong influence of local site effects on the amplitude of the seismic signal. Without station correction, the overall standard deviation of the magnitude residuals using our magnitude relation is 0.32, while using station corrections, the standard deviation 0.17. A comparison of our local magnitude relation with published mblg relationships for Southern Africa shows that these two magnitude scales are almost equivalent. The Ml equation for Central Southern Africa derived in this study has good correlation with mb values reported by the ISC and NEIC and attains systematically lower magnitude values than the South African relationship.


Introduction
Local magnitude scale (M l ) is one of the commonly used scales to quantify relative size of an earthquake. It was defined by Richter (1935) and Richter (1958) as the logarithm of the maximum zero to peak amplitude measured on a Wood-Anderson (WA) instrument with amplification of 2800 at a natural period of 0.8 s (Anderson and Wood 1925). Some studies (Uhrhammer and Collins 1990;Uhrhammer et al. 1996), as does the IASPEI standard formula (Bormann and Dewey 2012), have shown that the effective amplification of the typical WA seismograph is around 2080 times. The main aim of this study is to propose a relationship that estimates M l from the inversion of amplitudes of earthquakes recorded by broadband seismic stations located in Botswana, Zimbabwe, and South Africa. Results from this study will improve the calculation of earthquake magnitudes in the region. The study area is characterized by a moderate level of seismicity. Although, the area experiences, lowmagnitude seismicity, one large shallow earthquake occurred on the 3rd of April 2017 with a magnitude M6.5 (Asefa and Ayele 2021; Gardonio et al. 2018;Midzi et al. 2018). This earthquake was followed by aftershocks which lasted for several months (Mulabisana et al. 2021;Olebetse et al. 2020;Midzi et al. 2018). Monitoring of earthquakes in this area is carried out by the Botswana Geoscience Institute (BGI, formally the Botswana Geological Survey), Goetz Observatory (Department of Meteorological Services) in Zimbabwe, and also by the Council for Geoscience (CGS), South Africa. Seismic monitoring in Zimbabwe started in 1959 through the installation of the first seismic station in Bulawayo. The need to establish a seismic network in Zimbabwe was necessitated by the building of the world's largest man-made lake in Kariba. The current network is made up of four broadband seismic stations. One of these stations (BLWY) is operated by Africa Array (AF). In addition, MATP seismic station was installed in 2004 as part of the International Monitoring System (IMS) network. In 2014, Goetz Observatory installed 2 more seismic stations at Chipinge (CHIPN) and Karoi (KRI) to add to the existing stations.
Since April 1993 a global telemetered seismograph station (LBTB), located in Lobatse, has been operating in Botswana. From November 2013 to February 2018, a temporary array of seismic stations (Network of Autonomously Recording Seismographs (NARS), Botswana network) was deployed in Botswana to address the questions about the crustal and upper mantle structure. The NARS consisted of 21 broadband seismic stations distributed over the whole country. In February 2018, following the 3rd April 2017 M6.5 earthquake, the BGI converted the 21 temporary NARS stations to become part of the permanent network, which is now used to monitor earthquakes in and around Botswana. Data that is recorded from these stations is transmitted in real time to the BGI for analysis. Between 1997 and 1999, the Southern African Seismic Experiment (SASE) was set up and broadband seismic stations were installed in Botswana, Zimbabwe, and South Africa as part of this experiment. It was part of the Kaapvaal project (Carlson et al. 1996;Carlson et al. 2000) with the responsibility to investigate the seismic structure of the region. BGI, Goetz Observatory, and CGS use SEISAN (Havskov et al. 2020;Ottemöller et al. 2021) for manual analysis of local and regional earthquakes. BGI does not have its own magnitude relation and currently BGI is using the South African M l scale (Saunders et al. 2013). This is assumed to be reasonable since the two areas are close and tectonically similar. The geology of South Africa and Zimbabwe is similar; they both consist of Archean cratons which are almost of similar age and the cratons are both surrounded by mobile belts (Eriksson et al. 2011;Fouch et al. 2004;Jelsma and Dirks 2002). Similary, Goetz Observatory is using the California scale (Hutton and Boore 1987) which is not suitable for this area.
During the 1963-1991 period, the m blg magnitude scale was developed by Henderson (1974) and used for magnitude estimation in the southern African region by the Goetz Observatory. The Henderson (1974) relation (1) assumes shallow focal depth between 0 and 10 km.
in which A (max) is the maximum ground motion amplitude of the Lg phase in μm, T is the period in seconds, and the β[ ] is the distance normalizing term. Further details of this relation are found in Chow et al. (1980) and Hlatywayo (2001). The scatter in the derivation of β[ ] in Henderson (1974) was large and in order to improve the magnitude values of earthquakes in Southern Africa, Chow et al. (1980) developed an m blg magnitude scale using a large dataset and obtained the following distance normalizing term β( ) = 2.66 log 10 + 2.61 (2) for where A (max) is the ground motion amplitude of the Lg phase in nanometres (nm) and is the epicentral distance.

Data and processing
Three component digital recordings from 61 seismological stations from Botswana, Zimbabwe, and South Africa were used to carry out this study. These recordings correspond to shallow earthquakes of focal depth less than 35 km. Of the 61 stations used in this study, 20 are part of the SASE broadband seismic stations. Furthermore, we used seismic stations from NARS, BX-Botswana network and two Zimbabwean permanent stations. These stations continuously recorded earthquake data in that time period. Figure 1 shows a map of stations as well as epicentral locations of the 621 earthquakes whose parameters were used in this study.  (Saunders et al. 2013) The reviewed standard earthquake Bulletin that is routinely released by the International Seismological Centre (ISC) (publicly available at http://www.isc.ac.uk) was used to select 621 events for which matching waveform data were obtained from the Incorporated Research Institutions for Seismology (IRIS). Figure 2 shows the distribution of earthquakes used with respect to hypocentral distance. The analysis was restricted to earthquakes recorded in the range 0-1000 km from the hypocentre. Magnitudes in Fig. 2 were calculated using the relationship obtained in this study. Figure 2 also shows that most earthquakes used in this study are within the M l range of 2-7.
To generate an M l scale for Central Southern Africa, a Python-based algorithm was developed to automatically measure half peak-to-peak displacement amplitudes using the compiled waveform data. We took advantage of the software Obspy (Beyreuther et al. 2010), a Python-based package with several modules designed to help users in the workflow of downloading, inspecting and processing of seismic waveforms. All waveform traces were demeaned, filtered between 1 and 10 Hz, and transformed to displacement through deconvolution of the instrument response. Given the location of an event as provided in the ISC Bulletin, we used the AK135 velocity model (Kennett et al. 1995), to predict local P and S wave arrivals. We then applied a standard signal-to-noise ratio (SNR) detector around the predicted arrivals to check the quality of the waveforms (example of a processed waveform is shown in Fig. 3).
We performed visual examination on approximately 100 waveforms using different SNR's to identify the best minimum threshold for our study. We define the threshold to be 2.0; thus, waveforms with SNR ≤ 2.0 were rejected. For each station and for each earthquake, maximum trace amplitudes were measured on both N-S and E-W components. We returned the higher measured maximum amplitude between the 2 components as the amplitude of that station. The automatically measured amplitudes at most seismic stations compared favourably to amplitudes that are stored at the ISC Bulletin (e.g. Fig. S1). Our approach for amplitude picking described above is summarized in the flowchart in Fig. S2. The number Fig. 3 Regional earthquake seismograph 1997/08/01 4.7 m b(isc) for an earthquake located in South Africa that was recorded by station SA80 in Zimbabwe at 8.97 • from the epicentre. The top trace is the original record whilst the bottom trace is the simulated WA record, and p onset and s onset are the theoretical P wave and S wave arrival times respectively of observations made at each seismic station is shown in Fig. 4a and amplitudes picked are plotted against hypocentral distance as shown in Fig. 4b.

Methodology
A standard formulation of Richter's M l is written as (Richter 1935;1958;Hutton and Boore 1987): where A is half of the peak-to-peak amplitude (mm) of a horizontal component on a standard WA seismometer, − log A 0 is the distance normalizing term that reflects the overall attenuation attributes in the region of interest, and S is the station correction defined relative to a reference site condition. In order to maintain (Richter 1935) original definition of M l , − log A 0 is defined such that 1 mm of amplitude on a WA instrument located at a reference site 100 km away from an event would register as a magnitude 3 event.
The distance correction term is given as: The coefficient a depends on geometrical spreading, b on attenuation in the region of interest, combining Eqs. 3 and 4, the observed amplitude A ij in nm for each station for each event can be written as: Equation 5 is an overdertemined set of equation where i and j are earthquake and station index respectively. r ij is the hypocentral distance in km, S j is the station correction to the magnitude for the j th observation site, and C is a station correction term which gives the same magnitude for the same amplitude at the reference distance as in Southern California (Richter 1935). Equation 5 simply gives a system of linear equations which can be fed into a standard matrix formation in the form of: Taking p as the earthquake number and q as the station component number, Eq. 5 can be rewritten in matrix form as: − log( p(q−1) /100) −(r p(q−1) − 100) 0 · · · 0 1 0 · · · 0 −1 0 − log(r pq /100) −(r pq − 100) 0 · · · 0 1 0 · · · 0 0 −1 In Eq. 7, y ij = log(A ij ) + 3. The vectors of unknowns (x) from Eqs. 6 and 7 can be obtained by inversion of A using Eq. 8, under the constraint that the sum of station corrections is zero. The generalised inverse equation with at least N e + N s + 2 unknowns, represents a linear inversion problem that need to be solved using least-squares. N e and N s represent the number of events and stations respectively.
The linear system has a total of 684 parameters and 5306 data. The solution will be determined by singular value decomposition (SVD) method to estimate a, b, M l , and S j as implemented in the MAG2 program which is part of the SEISAN software package (Ottemöller and Sargeant 2013;Havskov et al. 2020;Ottemöller et al. 2021). We set the constant C to anchor the scale to the (Richter 1935) scale at 100 km. During inversion, an attempt to use a distancedependant a for up to three segments was attempted. However, this did not show any improvement over linear-segment hence the results are not discussed in this paper. Table 1 represents the input parameters that were used with the MAG2 code.
We calculated magnitude residuals between magnitude assigned to a single station and the median magnitude earthquake. The M l values obtained using the new relation derived in this study were then compared to those obtained for the same events using the Saunders et al. (2013) relation for South Africa, Hutton and Boore (1987) for California, and Ottemöller and Sargeant (2013) scale for the UK. Furthermore, our M l values were also compared to m blg values which we calculated based on the Henderson (1974) and Chow et al. (1980) relations. To evaluate whether the new scale is an improvement over any alternative relation, we look for a reduction in standard deviation (9) of all magnitude residuals computed as in which σ is the standard deviation of all the observations, N is the number of all observations, x i is the value of each observation, and μ is the mean of all the observations.

Results and discussion
The coefficients (a, b, and C) were determined by applying the SVD using the Numerical Recipe routines (Press et al. 1996) to invert 5306 amplitude value observations. The corresponding epicentres-station path covers most of the studied areas in Botswana, South Africa, and Zimbabwe (Fig. 5).
After inverting our observations, we obtained the following scale: M l = log 10 A+0.80×log 10 (R)+0.00086×R−1.37±S (10) in which A is ground displacement amplitude determined from instrument corrected synthetic WA seismograms in nanometres, R is the hypocentral distance (km), and S is the station correction. The estimated value for the geometrical spreading coefficient, a, is 0.80 ± 0.24, and the anelastic attenuation coefficient, b, is 0.00086 ± 0.00025 and C is −1.37. The obtained M l relation is valid for distances up to 1000 km. Figure 6 shows the comparison of our calibration curves against those for Southern California (Hutton and Boore 1987), South Africa (Saunders et al. 2013), and UK (Ottemöller and Sargeant 2013). It can be seen that for distances less than 80 km, the attenuation rate for Central Southern Africa is slightly larger than the (Saunders et al. 2013) relation, the (Hutton and Boore 1987) relation and the relation by (Ottemöller and Sargeant 2013). However, for all distances greater than 80 km, the attenuation obtained in this study is significantly lower than that of Southern California and the UK. This is not surprising since the geology of our study is dominated by old cratons which promote low seismic attenuation.
The station correction factor for a particular station represents the local site conditions (Richter 1958) and can be used to compensate for heterogeneous velocity structure near affected stations (Douglas 1967;Pujol 1988). High station correction values indicate site effects related to geological conditions or the presence of seismic noise have a strong influence on amplitude. Figure 7 is a representation of mean station residual variation across different networks within the study area. As we have explained before, BX network is a continuation of the NARS network. The only difference between the two networks is that during the transition, the sample rate was changed from 20 Hz (NARS network) to 40 Hz (BX network).
For the NARS and BX networks, stations are arranged such that the station below is a cosite station to the station above it. As an example, station NE201 is a cosite station to SKOMA and station NE213 is now station PHPEN. From Fig. 7, it can be seen that the two networks have similar station corrections. Also shown in Fig. 7 are the mean station residuals for cosite station BLWY an Africa Array (AF) station and SA72 from the SASE network. Station MATP is located at approximately 20 km south of station SA72/BLWY. As seen from Fig. 7, station SA72 has a low station correction and standard deviation. This is Fig. 7 Spatial variation of average station residuals (station magnitude-median magnitude) deduced in this study. Negative and positive average are shown by green and red colour, respectively. Text near the vertical bars denote seismic stations names. White polygons indicate the boundaries of the cratons inside our study area not surprising since it can be seen that almost all stations from the SASE network have consistently low station residuals as compared to the NARS and BX networks. Furthermore, the station correction differences between SASE stations and BX/NARS network stations is clearly visible at some nearby sites such as NE220 and SA68 as well as NE217 and SA71. We are not sure as to why the SASE network has low station corrections, but we can only speculate that this may be due to good metadata stored at IRIS. To sum up, station corrections calculated in this study range from −0.44 to +0.31 units of magnitude, which suggests significant differences in site effect.
We compared our newly developed M l , with magnitudes from other agencies. To do this, we used general orthogonal regression procedure (GOR) which accounts for errors in both variables (Gardner and Knopoff 1974;Pujol 2016). Compared to other regression techniques such as the standard linear regression (SLR), which is widely used, GOR is regarded as the best regression analysis procedure (Castellaro and Bormann 2007;Das et al. 2018). The ODRPACK which is a software package for weighting orthogonal distance from a set of observations was used for this task. The ODRPACK uses a FORTRAN-77 library that performs ODR with possibly non-linear fitting functions. The fitting functions are provided by Python functions operating on NumPy arrays. Figure 8 shows comparison of our M l with m blg derived from Chow et al. (1980), m blg after Henderson (1974), m b(I SC) , m b (NEI C) , and M l reported by South Africa M l(P RE) . To compare our relation with m blg , we extracted Lg maximum amplitudes from the ISC Bulletin. These Lg amplitudes were then used to compute the m blg magnitudes using relations after Chow et al. (1980) and Henderson (1974). For m b(I SC) , m b(NEI C) and M l(P RE) , we downloaded network magnitudes from the ISC Bulletin (Storchak et al. 2017; and then compared the network magnitude with event magnitudes calculated in this study. For our M l and m blg (Chow et al. 1980), we obtained the following GOR relationship: m blg (Chow) for m blg ≤ 5.9, n =100, Similarly, there is a good correspondence between the Henderson (1974), m blg and M l from this study. We derived the following relation   (H enderson) M l 0.91 ± 0.03 −0.11 ± 0.12 0.03 m blg (Chow) M l 0.87 ± 0.03 −0.26 ± 0.11 0.03 M (lpre) M l 1.01 ± 0.02 −0.62 ± 0.08 0.01 The model follows the relation y = α+βχ, where β is the slope of the line of best fit, α is the intercept, and σ 2 is the residual variance for m (blg) ≤ 5.9, n =100 For body wave m (bI SC) , the relationship follows: for m b(I SC) ≤ 6.7, n =90, Similarly, the regression between our M l and m b(NEI C) yielded the following relation.
for m b(NEI C) ≤ 6.8, n =84 For all GOR relationships discussed above, we obtained low residual variance between 0.01 and 0.04 (shown on Table 2).
Magnitude residuals, defined as the difference between the magnitude calculated at individual stations and the event magnitude for all the stations for a particular seismic event, are shown in Fig. 9 and Table S1. The magnitude residuals obtained using the relation by Hutton and Boore (1987) with station corrections is shown in Fig. 9a. Similarly, Fig. 9b shows results obtained using the relation by Saunders et al. (2013) with station corrections. Magnitude residuals obtained using the relation derived in this study without station corrections is shown in Fig. 9c and with station corrections in Fig. 9d. From Fig. 9, it is observed that there is less scattering in residuals for the newly developed magnitude scale with station corrections as compared to the residuals obtained using the other relations as well as the relation from this study without corrections. The overall, standard deviation of the magnitude residuals was computed to check the effect of the new magnitude relation on magnitude estimates. Summarised in Table 3 is the comparison of standard deviation of the magnitude residuals, and it can be seen that the standard deviation σ of magnitude residuals from our study without stations correction is 0.32. Similarly, with station correction, we obtained a standard deviation of 0.17. These results show that the addition of station corrections reduces the standard deviation and this reduction reduces station residuals close to zero. The standard deviation in residuals for Saunders et al. (2013) using station corrections is also quite low with values of 0.16. The average magnitude residual for different stations are shown in Fig. 9e for bin distances of 100 km, along with error bars. The residuals in Fig. 9e are small, which shows that the newly developed magnitude scale compensates attenuation effects for all distance ranges. Even though our new M l scale with station corrections and the Saunders et al. (2013) M l scale with station corrections produce almost identical residuals, a closer look at the number of observations (bar graph on right hand side in Fig. 9) indicates that the average residual of our results with station correction is close to zero compared to the Saunders et al. (2013) relation. The geometrical spreading parameter and the distance attenuation parameter values are slightly different.

Conclusion
A local magnitude scale was derived for Central Southern Africa using 5306 traces of maximum amplitude observations. The observations were measured on the horizontal components of observed WA seismograms at distances ranging from 0 to 1000 km. The data used was obtained from records of 621 shallow earthquakes recorded by 61 broadband seismic stations. The scale derived was based on the (Richter 1935) definition of M l . As part of the inversion process, station corrections were determined for each of the stations that contributed data. Values determined range from −0.44 to +0.31. Magnitude values obtained using the new Central Southern Africa M l scale were compared to those published by NEIC, PRE and ISC. The comparison shows that our solutions are stable and can replace the current M l relation being used in the region. Results obtained were also analysed by studying the trend of the magnitude residuals. Overall, the standard deviation of the magnitude residuals without using station corrections is 0.32 Table 3 Comparison of the geometrical spreading parameter a, distance attenuation parameter b, and standard deviation σ in magnitude residuals of the Saunders et al. (2013)  while those obtained in this study using station corrections is 0.17. The reduction in standard deviation shows that using the new relation obtained in this study brings the average residual close to zero.

Data and resources
The bulletin data that was used in this study were downloaded from the ISC Bulletin (ISC 2022). Waveform data were obtained from the Incorporated Research Institutions for Seismology (IRIS). The python package Obspy (Beyreuther et al. 2010) was used to obtain some of the waveform data. Most of the plots were generated using the Generic Mapping Tool (GMT) of Wessel et al. (2019). Figure 8 was generated using Matplotlib (Hunter 2007)