Re‐examination of the 1928 Parral, Mexico earthquake (M6.3) using a new multiplatform graphical vectorization and correction software for legacy seismic data

Modern analysis of earthquakes depends on digital time series; however, only about 30% of the timespan of recorded seismicity is available in digital format. During the first half of the 20th century, most of the earthquake ground motions in Mexico were recorded by Wiechert mechanical instruments on smoked paper. In this work, we developed and use Tiitba, a new portable multiplatform graphical user interface (GUI) open‐source software coded on python, to vectorize and correct old analogue seismograms. Using this software, we vectorized the 11/1/1928 Parral (M6.3) earthquake seismograms to obtain the constant time‐interval digitized time series for each component and station of the seismogram. Then, we obtained its source focal mechanism using a genetic algorithm methodology. Also, with an auxiliary Tiitba module, we constructed a SEISAN S‐file to relocate the hypocentre of this earthquake. The obtained relocation is about 125 km south of a recently recorded seismic swarm that occurred in 2013, which has a similar focal mechanism to the one that we obtained in this work.


| INTRODUCTION
Digitizing analogue seismograms has become a necessity to preserve historical earthquake information and to perform modern seismicity re-analyses with the recovered digital time series seismograms. Analogue seismograms represent about 70%-75% of the timespan of instrumental seismology (e.g. Wang et al., 2019;Wang, Jiang, et al., 2016). Thus, these seismograms are of immense scientific value.
During the first half of the 20th century, most of the earthquake ground motions were recorded by mechanical instruments on smoked paper, photographic film or ink on a whitepaper. After early 1970, standardizing efforts such as the microfilms distributed by the World-Wide Standardized Seismograph Network (WWSSN) and other initiatives around the world helped to give better access to seismological information and seismograms (e.g. Batlló et al., 2008).
In different published works on historical earthquake re-analyses, authors have used manual or semiautomatic vectorization schemes to obtain digital time series from analogue seismograms (e.g. Baştürk et al., 2016;Ichinose et al., 2003;Pino et al., 2004, among others). Different GUI-based software to obtain digital seismograms from raster images have been developed, for example, SeisDig by Bromirski et al., (2003) and DigitSeis by Bogiatzis and Ishii (2016), both developed under the MatLab platform. Also, the TESEO programme as part of the SISMOS project by Pintore et al., (2005) developed as a GIMP-GNU plugin, as well as other additional programmes, for example, Church et al., (2013), and , Wang, Jiang, et al. (2016). For the case of old or partially damaged smoked seismograms with very thin and discontinuous traces, usually the use of semiautomatic functions for vectorizing these records results in time-consuming procedures requiring additional tasks.
Most of the seismological stations of the Mexican network during the first half of the 20th century were equipped with Wiechert mechanical seismographs recording on smoked paper. The Wiechert seismograph records commonly contain remarkably misaligned time-marks among the traces, very thin and discontinuous traces, relatively damaged smoked papers, and different inhomogeneities in the smoking densities in the paper. Because of this, in order to facilitate the vectorization and correction of the old and relatively damaged smoked seismograms from this Mexican network, we developed Tiitba. Tiitba is a graphical user interface-based software developed ad hoc to vectorize and correct old smoked paper seismograms. Using this software, we vectorized and corrected the seismograms from the 11/1/1928 earthquake that occurred north of Parral city in Chihuahua, Mexico. This earthquake is one of the largest intraplate cortical earthquakes recorded in northwest Mexico in the last century. Parral earthquake has been previously analysed by Doser and Rodriguez (1993) where they computed the focal mechanism by inversion of teleseismic records, constraining the procedure by the first motion polarities of the data and by the characteristics of the faults in the region. The 1928 Parral earthquake epicentral zone is located south of a seismic swarm (Mw5.3) which occurred in 2013, showing evidence that this represents an active region of the Basin and Range province of Chihuahua with extensional deformation stresses (Suárez et al., 2016).
Using the vectorized and corrected time series, we recomputed the source parameters through a genetic algorithm procedure. Also, with an auxiliary Tiitba module, we construct a SEISAN S-file to relocate its epicentre.

USER INTERFACE (GUI) FOR HISTORICAL SEISMOGRAMS ANALYSIS, VECTORIZATION AND CORRECTION
Tiitba (the pre-Hispanic Mayan word for 'the earth shakes') is an interactive open-source GUI developed and completely coded on Python. This software can be executed in all Python and C+ compatible operating systems (Windows, iOS, Linux, iOS and Android). In Tiitba PyQt ® is used as the graphical link between the code and the operating system, employing the Open-source Computer Vision Library (OpenCV) for the high-resolution raster images visualization and manipulation; a detailed description of the libraries, coding and references for its use is included in Appendix S1. Tiitba GUI is structured in three main modules and two auxiliary functions modules, each of them with specific procedures to analyse, vectorize and correct the analogue seismograms.
In Figure 1, we show the GUI main window, which is divided into three sections: the first one for function buttons, the second one for the information panel and the third one for showing the graphical representation of Matplotlib canvas (Hunter, 2007).

| Digital image processing (DIP) module
Wiechert seismograph historical records commonly contain strongly misaligned time-marks among the traces, marginal handwritten notes, annotations relative to the corresponding seismic events and different inhomogeneities in the smoking densities in the paper. Thus, in most cases, digital image processing of the scanned seismograms is required before vectorizing.
This module involves the procedures for image trim, contrast increase and conversion into a binary black/ white image to highlight the seismogram trace from the background, either smoked (black) or white paper. After processing, the image disc size can be also reduced up to 80%.

| Seismograms vectorization (SV) module
There are at least four common general problems when vectorizing old analogue seismograms. The most common arises from the graphical interruptions of the seismogram traced lines, especially for large seismogram amplitudes. A second problem involves the variations in the smoke densities and inhomogeneities of smoked areas in analogue F I G U R E 1 Tiitba GUI (a) Main Window, functions modules, panel information and Matplotlib canvas for graphical representation. (b) Auxiliary modules window, the left tab is for phase picking function to create a SEISAN S-file, and the right tab is for improvements to time series data functions seismograms due to the inherent irregularities during the smoking process. A third problem derives from the trace amplitude saturation produced by the restricted dynamic range of the mechanical seismographs. Finally, the continuous intersection of the traced lines with the previous and subsequent traces of the seismogram, especially for large regional earthquakes recordings, usually masks and even conceals important amounts of information. Most short-period ground motion records produced by large earthquakes from Wiechert seismographs have crossing traces, totally clipped amplitudes, and/or records are incomplete due to the dislodgement of the stylet. All these common problems strongly complicate following the trace for correct vectorization.
The SV-module simplifies the panning and zooming over the raster image to follow in a chronological sense the seismogram trace and vectorize the record for as long as possible, despite the imperfection issues.
This module has a raster image size function, which computes the analogue paper size with the density value on the metadata raster image in pixels per inch (PPI). This module includes two different functions for the scaling of the output time series: one by picking at least three continuous timemarks on the seismogram image and one by inputting the time and amplitude values at opposite corners. Once these values are set, the vectorize function is employed to convert seismograms from analogue to digital format. The module also included the plot and save data functions.

| Basic corrections (BC) module
Raw data from the manually vectorized seismogram do not have precise evenly spaced time samples. These data also contain trends and curvatures due to the mechanical instrument drawing device. Therefore, the main corrections for these data depend on the type of seismograph and some specific technical characteristics of each seismogram, which can be grouped into the geometrical corrections, timing corrections and instrument corrections (e.g. Batlló et al., 2008).
In this module, it is possible to apply basic corrections such as: time series detrending, time resampling for a constant (uniform) time sampling interval, curvature effect removal and additional corrections for the instrumental response (Wiechert mechanical seismographs).
Tiitba modules, procedures and functions are detailed in Appendix S1.

| Auxiliary modules
Technological and scientific developments allow performing a modern re-analysis of historical earthquakes data, either from vectorized time series ASCII files or from raster images (Figure 1b). Taking advantage of these advances, we develop two auxiliary modules for historical earthquakes relocation and time series improvement.
The module for earthquake relocation allows picking of the P and S phases directly on analogue seismogram images to create a SEISAN S-file (Havskov et al., 2020) for hypocentral relocation. As in SEISAN, in this module, the user can add the weight and select the arrival type (either impulsive or emergent) for each phase.
The module for time series improvement consists of an ad hoc algorithm of weighted least squares to recover clipped amplitudes. The algorithm is applied for each clipped pick.
In this auxiliary module, a function is also included for saving the data in at least two standard seismic formats either SAC or miniSEED by entering the date and time as well as the station and the instrument metadata.
For more details on the Tiitba GUI modules and functions, refer to Appendix S1 of this work.

EARTHQUAKE, CHIHUAHUA, MEXICO
The M6.3 Parral earthquake occurred on 1 November 1928 north of Parral city in Chihuahua north-western Mexico ( Figure 2). Doser and Rodriguez (1993) computed the focal mechanism of this earthquake using teleseismic P-wave data, constraining the solution with the fault parameters from different mapped faults in the epicentral area. They conclude that the rupture characteristics of this earthquake are similar to previous seismic events of comparable magnitude occurring within the western Cordillera of the USA (Doser & Rodriguez, 1993). More recently, a seismic swarm, 125 km north of the Parral earthquake, was well recorded in 2013, where the largest event (Mw5.3) occurred on September 21 in the same clustered epicentral zone. Suárez et al., (2016) located and analysed more than 200 events of this swarm and computed the focal mechanism of five of the largest events through a temporal seismic network. They concluded that most of the five focal mechanisms from the 2013 swarm are highly similar to those obtained by Doser and Rodriguez (1993) for the 1928 earthquake, suggesting an ongoing northeast-southwest extension of the Chihuahua Basin and Range region.

| Legacy data and metadata
The Parral 1928 earthquake was well recorded by the National Seismological Service of Mexico (SSN) seismographic stations at regional distances. Nowadays, the Mexican analogue seismographic collection is stored at the Joint Library of Earth Sciences of UNAM (BCCT), as part of the Sismoteca-SSN project to preserve, digitize and reuse this amount of legacy data (Pérez-Campos et al., 2020). Due to the availability of the Sismoteca-SSN seismographic collection, we have been able to get and use the original 1928 seismograms from the Mazatlan (MZX), Guadalajara (GUM), Merida (MER), Tacubaya (TAC) and Manzanillo (MNZ) stations. Each analogue smoked seismogram was scanned at 1,200 DPI and saved in JPG image format resulting, for an average 90 cm length seismogram, in a disc size for each image of about 40.0 Mb. We obtained the arrival times, instruments constants and earthquakes observations from the 1928 Mexican earthquake's bulletin (Servicio Geologico Mexicano, 1930). From the Caltech online archive for analogue seismograms (Brodsky & Kanamori, 2011), we also were able to get the scanned analogue seismograms from the Riverside (RVR), Pasadena (PAS) and Santa Barbara (SBC) stations. California's first seismographic stations mainly include Wood-Anderson seismographs, and from the Brodsky and Kanamori (2011) scanned images of the Caltech Archive report, we obtained the seismographs metadata. Since all these three stations lie along a very similar back azimuth with respect to the Parral earthquake epicentre, we use only RVR station due to the highest amplitude and better signal-to-noise ratio relative to the other two.
In Table 1, we present the available data for location, magnitude, focal mechanism and other seismic information for the 11/1/1928 Parral earthquake. Additional metadata, data and seismic information can also be found in the literature (e.g. km McComb & West, 1931) and on online bulletins (e.g. International Seismological Centre, 2021, online Bulletin).

| Data processing
For the SSN stations, we only process the data from the closest stations to the epicentre with the highest amplitude and the best signal-to-noise ratio in the first arrivals.
For the Caltech stations, we only process the data from RVR station. Using the Tiitba DIP module, we perform the digital image processing of the seismograms from the MZX, GUM and RVR stations, to highlight the seismogram trace from the background for vectorization. We used the SV-module to vectorize these seismogram images. First, the amplitude and time output scale were obtained using the continuous time-marks options. Then, we redraw the seismogram traces and finally saved the sampled traces as a two-column ASCII file, where the first column corresponds to the uneven time samples in seconds, and the second column corresponds to the amplitudes in mm (based on the original analogue paper seismogram size).
To obtain the constant time-interval velocity time series, we used the BC-module. First, we detrended and corrected the seismograms to their baseline. Then, we corrected the SSN seismograms for the curvature effect produced by the recording drawing stylet on the Wiechert seismographs, fixed on one side to the seismograph and damping system. Tiitba curvature correction is based on the Grabrovec and Allegretti (1994) procedure (see Appendix S1), which requires two essential inputs from the metadata archives: first, the stylet length, which has been directly measured (physically) on the original Wiechert seismograph at the TAC station currently belonging to the Museum of Geophysics of UNAM; and second, the recording drum speed, which is computed using the output scale function. Next, we resampled the seismogram traces at 20 Hz sampling rate. Finally, we performed the instrumental correction. For the SSN records, we used the MC-module which requires different seismographic metadata inputs: the magnification (V ), the damping rate constant (ε) and the seismograph natural period (T o ), which we obtained from the 1928 Mexican earthquake bulletin by Instituto de Geologia (1930). For the Caltech Archive seismograms, we computed the Wood-Anderson instrumental response using an ad hoc Python code, the pendulum's natural period T and damping constant ε were obtained from the Scans of the Caltech Archive report (Brodsky &Kanamori, 2011) andH. Kanamori (personal communication, 2021). Table 2 lists the correction parameters and constants used in this work. Figure 3 shows the final processed seismogram images and the velocity time series for all the SSN and Caltech seismograms recordings used for this work.

| Hypocentral relocation
One of the most used and tested software for earthquake location and relocations is SEISAN (Havskov et al., 2020). Tiitba GUI has an additional module that allows picking the P and S arrival time directly from the analogue seismogram images and then saves the information in an S-file format.
In the SSN seismogram raster images of the 1928 Parral earthquake, the reference time-marks were defined on the nearest minute relative to the origin time reported by the SSN. Due to the unsynchronized clocks for the mechanical seismic stations in 1928, some of the absolute travel times calculated from the recorded times at the different stations do not have real physical meaning. Therefore, for the 'P' phase the pick is assigned in this case, a weight in the S-File of value 9, which means that in SEISAN only the S-P time is used in computing locations. We pick 'P' phase arrival on the vertical seismogram and 'S' phase arrival on either of the two horizontal seismograms (Figure 4).
To relocate the earthquake hypocentre, we used SEISAN and an updated velocity model for Mexico used by the SSN (Table 3). We used the geographic locations of the 1928 seismic SSN national network. After the relocation, the hypocentral location obtained in this work is Lat = 27.333°N, Lon = 105.464°W, Depth = 10 km and an RMS = 0.47 error.

| Focal mechanism
Using the new location result, and the vectorized and corrected P-wave time series obtained here, we performed a global search for the focal mechanism through the Pygad genetic algorithm (GA) code with 4 degrees of freedom, entirely coded in Python (Gad, 2021). For this work, we used the GA algorithm for the global solution, the complete angle range for each focal mechanism parameter (for the strike angle φ = [0°, 359°; Δφ = 1.0°], for the dip angle δ = [0°, 90°; Δδ = 1.0°], for the rake angle λ = [−180°, 180°; Δλ = 1.0°] and depth D = [0.0 km, 12.0 km; ΔD = 2.0 km]) through a fitness function that compares the observed data with the synthetic data. The fitness function used here for the global search in the GA was the semblance function shown by Equation (1): where S(m) is the semblance between observations and synthetics for each component, cross d s d o is the crosscorrelation between synthetics and observations, autoc d s is the autocorrelation of the synthetics, and autoc(d o ) is the autocorrelation of the observations. Then, as S(m) becomes closer to zero, the similarity increases.
For the direct model, we used the Discrete Wavenumber Method programme for synthetic seismograms AXITRA (Cotton & Coutant, 1997), and the 1D velocity model for computations was obtained for each source station path, based on the 3D S-wave velocity model for Mexico and South US by Spica et al. (2016) ( Table 4).
To test the accuracy of our method, we first computed the focal mechanism for the 21 September 2013 (M w 5.3) earthquake and compared it with the reported one by different authors. For this, we used the P-wave phase data recorded at the Topolobampo (TSIG), Hermosillo (HSIG) and Mazatlan (MAIG) broadband SSN stations (Figure 2, red triangles) and the seismic moment reported by Franco et al. (2020) from the Moment Tensor Catalogue for Mexican earthquakes (MTCME).
Our inversion result using the GA global search method for the 2013 earthquake is = 185 • , = 87 • , = − 81 • and depth D = 5.0 km), with a mean semblance after inversion accounting for all the three orthogonal components of S(m) = 0.219 ( Figure 5). This result agrees well with the mechanism reported by the MTCME ( = 196 • , = 89 • , = − 77 • and D = 2.0 km), using the same type of station's data (regional broadband recordings; Figure 6). However, this result becomes not so similar compared to the one reported by the Global CMT catalogue (Dziewonski et al., 1981;Ekström et al., 2012) ( = 22 • , = 59 • , = − 29 • and D = 13.2 km); it is to be noted that for this result they used only teleseismic long  period waveforms. To quantify the difference between our focal mechanism result and the one reported by the MTCME and Global CMT, we computed the Kagan angle K, being this the minimum 3D angle required to rotate the principal axes of one of the moment tensor onto the other (Kagan, 2007): K = 0 • means that both nodal planes match exactly, otherwise K tends to 120° at the maximum difference. The Kagan angle between our solution and the MTCME solution is K = 11.43° and between our solution and the Global CMT is K = 75. 52 • .
F I G U R E 5 Focal mechanism result, by the genetic algorithm method for the 09/21/2013 Earthquake M = 5.3, Chihuahua northeast Mexico. (a) Shows a box plot for the source parameters misfit, where ϕ is the strike in a range between 1° and 360°, δ is the dip in a range between 0° and 90°, λ is the rake or slip in a range between −180° and 180°, and D is the depth in a range between 1 and 11 km. (b) Shows the semblance mean decay over the genetic algorithm generations. (c) Shows the fit between the data (solid black line) and the synthetics (dashed blue line) and the focal mechanism solution. The semblance mean is 0.219, and the Kagan angle K is 11.43° Figure 6 shows the result of the global search for the focal mechanism of 2013, M w 5.3 earthquake, compared with the published result by (Franco et al., 2020) and the Global CMT catalogue (Dziewonski et al., 1981;Ekström et al., 2012).
For the computation of the focal mechanism of the 1928 Parral earthquake, we used the data obtained and processed in this work from MZX, GUM and RVR recorded on smoked paper. In this case, we used only the initial 20 s of the P-wave phase from each station because this phase, in general, maintains its complete amplitude and generally is unclipped during the initial ground motion. We execute the GA code by also changing the seismic moment estimates within the interval reported by Doser and Rodriguez (1993) 3.2 ± 1.0 × 10 25 dyne-cm.
Results from this inversion are shown in  Figure 8 shows the comparison between the focal mechanism results. Doser and Rodriguez (1993) used teleseismic data and a restricted range for the plane parameters. In this work, we used regional records and F I G U R E 6 Comparison of the focal mechanisms for the 09/21/2013 Chihuahua earthquake M = 5.3, Northeast, Mexico, for the solution reported by Franco et al., (2020) in the MTCME, the one reported in the Global CMT catalogue, and the one obtained in this work. Coloured symbols indicate the earthquake epicentres reported by different published works and seismological agencies. The yellow star represents Parral town considered a complete global solution search. Although the Kagan angle suggests that the nodal planes between solutions are not close to each other, the parameters misfit in our solution indicate a better agreement when a global search range is used (Figure 7-1).
Our result for the 1928 earthquake shows a normal faulting mechanism with a strike-slip component. Both nodal planes should be considered as the fault rupture planes; however, the one striking north-south, which agrees better with the strike angle of the mapped faults in this area, F I G U R E 7 Focal mechanism result, by the genetic algorithm method for the Parral earthquake 11/01/1928 Earthquake M = 6.3, Chihuahua northeast Mexico. (a) Shows a box plot for the source parameters misfit, where ϕ is the strike in a range between 1° and 360°, δ is the dip in a range between 0° and 90°, λ is the rake or slip in a range between −180° and 180°, and D is the depth in a range between 2 and 12 km. (b) Shows the semblance mean decay over the genetic algorithm generations. (c) Shows the fit between the data (solid black line) and the synthetics (dashed blue line), and the focal mechanism solution. The semblance mean is 0.139, and the Kagan angle K is 52.61° appears more likely to be the rupture plane. Also, this mechanism is very similar to the focal mechanism of the 2013 earthquake ~125 km north from our relocation.

DATA IN MEXICO
For more than half a century, the Mexican seismic data were limited to no more than 9 seismic stations. In most cases, all the recovered and available seismic data have been classified and appropriately stored at the BCCT. The SISMOMex project (Pérez-Campos et al., 2020) has successfully preserved and stored a large number of analogue seismic records from the early 20th century to 2015. For the large earthquakes that occurred during the first half of the 20th century, the Mexican seismic records were limited to nine stations at the best case. When aiming to reuse this very old data, and towards a re-analysis of these large earthquakes, we face limitations regarding the azimuthal coverage on a regional scale. Nevertheless, the legacy data are fortunately well preserved and now available for its reuse. Nowadays to perform a re-analysis on very old earthquakes, it is necessary to use all legacy data available, taking advantage of the current state of the art methodologies to expand our knowledge on the local and regional seismicity and its tectonic implications. The SSN and BCCT are now preserving Mexico's legacy data, as F I G U R E 8 Comparison of the focal mechanisms for the Parral earthquake M = 6.3 11/01/1928 Chihuahua, Northeast Mexico, for the solution reported by Doser and Rodriguez, (1993), and the one obtained in this work. Coloured symbols indicate the earthquake epicentres reported by different published works and seismological agencies. The yellow star represents Parral town well as bulletins, and linking this information to the SSN repository. Seismograms are being carefully scanned and metadata linked to the image, and digital information has begun to be incorporated (Pérez-Campos et al., 2020).

DISCUSSION
In this paper, we use Tiitba, a new portable multiplatform software written in Python, for vectorizing the analogue seismograms of the Parral 1928 (M6.3), Mexico, earthquake. Detailed information about the Tiitba-GUI coding, specifications and use can be found in Appendix S1.
Using Tiitba, we obtained the vectorized time series for regional legacy seismic seismograms recorded on smoked paper and ink on paper for the 1 November 1928 Parral earthquake M = 6.3. Then, we corrected these time series by detrending, correcting for curvature and instrumental response and performed the focal mechanism inversion by global search using a genetic algorithm method, employing a specific semblance fitness function. The results show a normal focal mechanism with a strike-slip component (ϕ = 180°, δ = 20°, λ = −50°). The mechanism is highly similar to the mainshock of 21 September 2013 (M W = 5.3) swarm, which has been one of the largest, well-recorded swarm in the Chihuahua Basin and Range region; this earthquake occurred approximately 125 km north of the Parral earthquake epicentre and may be considered an extension in this region of the Basin and Range (Suárez et al., 2016).
Using an auxiliary module of Tiitba, we obtained the P and S phases' arrival times from the 1928 earthquake at the different regional stations, directly on the analogue seismograms, and generate a SEISAN S-file to relocate the earthquake. The relocation of this hypocentre shows that the earthquake occurred 26 km north from the locations reported by Sawires et al. (2019) and the ISC-online bulletin (2021). Our results show an acceptable RMS location error, considering the complexities involved in obtaining the seismological information from the early 20th-century analogue seismograms, especially those related to the absolute times of the P and S phase arrivals.