On the recovery and analysis of historical seismograms

The analysis of historical seismograms has proven to be a fundamental tool to help with the definition of the seismic risk in specific regions. Indeed, modern quantitative reappraisal of relevant earthquakes that occurred before the 1960’s; i.e., prior to both the developments of modern recording instruments and the theoretical progress, has been essential for the assessment of the seismic potential of a source area. However, due to the characteristics (transducing and recording) of the old analog seismographs, the data available are affected by intrinsic uncertainties, and errors can be introduced during the processing of waveform digitization. These drawbacks can seriously influence the quality and reliability of an investigation. In general, no standard technique can be applied when dealing with historical seismograms. Thus, specific tests and cross-checks have to be designed to estimate the limits of each specific analysis. Here, we aim to provide an overview of the whole procedure while focusing on the most crucial steps, from the seismogram recovery to the application of modern techniques for the retrieval of the seismic source information. We also suggest possible checks for the robustness of the data and for the available instrument characteristics, with a description of the effects of various uncertainties on the results that can be obtained. We thus provide useful indications for the analysis of historical seismograms, and also for the correct interpretation of the resulting characteristics of the seismic source.


Introduction
Major faults represent stable geological structures where huge amounts of strain energy can accumulate in response to stress produced by tectonic movements.This energy is then released in a few seconds, giving rise to major earthquakes, which is a process that occurs recursively, and often in a similar way on the same fault.As it takes a long time to accumulate large amounts of energy, the biggest earthquakes have long return periods, which range in many cases from several hundreds to thousands of years.Knowing the energy released by major earthquakes, and the slip kinematics, size of the ruptured surface, and distribution of displacement on the fracture is fundamental for the definition of the potential of the causative seismogenic faults, and thus for the assessment of the associated seismic hazard.
Currently, seismic source studies are mainly based on the analysis of the large amounts of digital data that are collected with modern instruments.These data allow for detailed reconstructions of a seismogenic processes.On the other hand, modern instrumental seismology was developed only about 50 years ago, so since the 1960's, when modern electromagnetic seismometers and, shortly later, digital recorders were designed and installed all over the world.With recurrence times of hundreds or thousands of years, detailed information on the characteristics of an earthquake source is only available for a relatively small number of significant events.On these grounds, investigations into historical earthquakes are fundamental, and, indeed, they have become progressively more important in recent years.
Overall, research into historical earthquakes can be grouped into three main distinct fields that are characterized by the differences in both the types of data and the methods of investigation: paleoseismology, macroseismicity (studies of intensity data) and analysis of coseismic geophysical data.This last approach is based on objective measurements of variables that relate to the whole space-time evolution of the fracturing processes, and it represents the most appropriate framework for analysis of the characteristics of a seismic source.
Usually, instrument data of historical earthquakes comes from geodetic measurements of static surface displacement or from recordings of transient ground motion.To quantitatively estimate the permanent surface deformation generated by a seismic event from geodetic data, a reference measure in the source area should to be available, which should be taken from before the earthquake happened [e.g., Bilham andEngland 2001, Nyst et al. 2006].Unfortunately, this is not often the case, and in particular for events older than a few decades.
Unlike static displacement, seismic waves generated by strong earthquakes produce effects that are detectable over the entire planet.Thus, considering that since the beginning of the 20 th century a relatively large number of seismic stations have been operating globally, historical seismographic data can potentially provide a source of precious data for investigations into major seismic events that occurred before the electronic instrument era.Indeed, depending on the quantity and geographical distribution of the recording stations at the time of a specific earthquake, there might be a sufficient number of seismograms to be used to reconstruct the associated process of rupture.
However, quantitative seismological analyses and modern perspectives of a seismic source that allow the description of a rock fracture in the Earth crust starting from the analysis of the recorded seismograms were only developed from the early 1960's.Thus, to gain more information on a source of an important earthquake that occurred from 1900, the original historical seismograms must be gathered and re-analyzed.Even though they are analogous to modern recordings, at least in principle, these seismograms have fundamental differences relative to digital waveforms, and special care is needed for their collection and analysis.
Several studies have been published that have described the general procedures for quantitative analysis of historical seismograms [e.g., Batlló et al. 2008, Pino 2011], which have also shown several case histories.In the present study, within a general overview on the recovery, processing, and study of historical seismograms, we focus on some of the possible difficulties, the results attainable, and the intrinsic limitations of the analysis.The aim is to provide both useful indications for these kinds of investigations, and valuable information for the interpretation of the results of the analysis of historical seismograms.

Retrieval of historical seismograms
The first step in the study of historical seismograms is the data retrieval.Since the early 20 th century, the total number of seismological observatories worldwide greatly exceeded 100 [Schweitzer and Lee 2003], and many sites have operated more than a single instrument.At the same time, the number of seismograph stations rapidly increased in the following years [Johnston 1996].Seismographs were based on either mechanical or electromagnetic equipment, and up to the last decades of the 20 th century, analog data were recorded on plain, smoked or photographic paper.In most cases, the paper sheets were changed daily, therefore a large amount of paper began to accumulate at observatories very rapidly, which became a problem due to both the need of suitably large spaces for storage and the perishable nature of the support.Over time, with decreasing interest, with wars, and with the closing (or relocation) of observatories, these original recordings have been lost and damaged.Moreover, at that time, the analyses were performed on the original recordings, and in particular, for the most significant earthquakes the seismograms were often loaned all over the world, and were not always returned to the original observatory.Consequently, many recordings were destroyed or lost.
Fortunately, in recent decades, the scientific community has become increasingly aware of the importance of preserving this huge amount of data, and several international initiatives have been proposed with this aim [Lee and Benson 2008].In particular, specific projects for the collection, recovery and storage of historical seismograms have been implemented [e.g., SeismoArchives, http://www.iris.edu/seismo/projects;SISMOS: Michelini et al. 2005;EuroSeismos: Ferrari and Pino 2003], which have generally made the scanned images available.
To be analyzed by means of modern techniques, historical seismograms need to be converted into waveforms that are analogous to present digital recordings.This process is generally accomplished through separate, successive steps, starting with the scanning of the original image.Thus, waveform digitization can be performed with automatic or semi-automatic procedures.Several studies have developed computer programs for the digitization of historical seismograms [e.g., Baskoutas et al. 2000, Bromirski and Chuang 2003, Pintore et al. 2005] and more codes are still being developed [S.R.Taylor and X. Yang: http://www.rdss.info/librarybox/mrr/MRR2010/PAPERS/08-07.PDF].
Digitization techniques based on the manual selection of the points are mainly characterized by distinct methods of interpolation between the single points, to obtain a waveform that can then be resampled at a constant time interval.A crude linear interpolation would give a polygonal line that would account for the main signal features, but would also strongly alter its frequency content.On the other hand, more sophisticated interpolation methods that are based on the global properties of the sample data (such as splines or Fourier-based techniques) or on the local characteristics of the signal (such as those proposed by Akima [1970] or Wiggings [1976]) would reproduce the signal more faithfully, although a careful check is still required to avoid undesired oscillations that are not present in the original recordings.As for automatic or semi-automatic techniques, the start and end points are usually selected, and thus the single points in-between are searched according to a chosen sampling rate.This process is generally based on color analysis of the images [e.g., Bromirski and Chuang 2003] or even making use of neural networks [e.g., Pintore et al. 2005].Overall, the definition of a best technique is not possible, as it will depend on both the features of each specific seismogram and the expectations of the operator.
Waveform digitization is a critical step in data processing, whether it is manual or automatic: signals generated by actual ground motion can be taken out, or conversely, inconsistent features can be introduced that were absent in the recording.Such drawbacks are more frequent when the original seismogram is not very clear or is partially missing (e.g., signal thickness is too large compared either to the wave amplitude or to the distance between successive oscillations; holes in the waveform due to time marks; very fast oscillations of the needle that prevented clear writing) (Figure 1).
Among those especially developed for seismic waveforms digitization, most codes also allow for the correction of possible distortion due the finite arm length and inclination (Figure 2), with the transformation of each generic plane curve obtained from the vectorization into a single valued time function for the amplitude.In general, these corrections are made with numerical algorithms that require the estimation and/or assumption of some parameters, thereby introducing elements of subjectivity.Figure 3 shows a comparison of two seismograms derived from the same original recording and processed by two different operators using the same procedures.The difference between the waveforms is evident and becomes macroscopic when a low-pass filter is applied to the seismograms.These differences can have important effects on the results.
Finally, when the original recordings are being transformed into digital waveforms, some caution must be exercised with the determination of the paper speed.This can be achieved by carefully checking the distance between the time marks, as historical seismographs can sometimes experience significant changes in paper speed over a few minutes.

The analysis of historical seismograms
Although there are the limitations and possible errors as described in the previous section, seismograms from this processing are analogous to the signals recorded by modern digital stations, in that they represent traces recorded by seismographs in response to ground motion at specific sites.Thus, modern techniques of waveform analysis and modeling can be applied to the study of the seismic source and to the structure of the Earth interior.
The shape and the amplitude of a seismogram is determined by the transduction characteristics of the instrument; i.e., the way in which the chain formed by the sensitivity to ground motion, the damping apparatus, and the writing system transforms the oscillations produced by the seismic waves.In modern seismometers based on electronic technology, the instrument response is determined by the characteristics of the circuit components that make up the various parts of the instrument, and these are well known and extremely stable, and can be verified with accuracy.On this basis, the actual ground displacement can be reconstructed with negligible error.The analytical form of a frequency response can be derived from the instrument characteristics.For example, in the case of an old mechanical instrument (such as a Wiechert seismograph), this depends on the natural period of the pendulum, the damping coefficient, and the static magnification (electromagnetic seismographs, like Galitzin seismographs, also have a galvanometer natural period).The instrument parameters are usually reported in the station logs or in the periodical station bulletins, and so this is why these records are extremely important: no waveform analysis is possible without knowing the characteristics of the seismograph.

ON THE RECOVERY AND ANALYSIS OF HISTORICAL SEISMOGRAMS
However, even when these bulletins are available, in some circumstances problems can arise from possible instrument instability.Indeed, especially for seismographs used during the last decades of the 19 th and early 20 th centuries, the responses to ground motion were very unsteady, and usually the parameters were not tested very frequently.At Potsdam, for instance, the natural period of the horizontal Wiechert seismograph that corresponded to 14 s in November 1905 for both components, turned out to be 19.8 s and 12.5 s for the N and E components, respectively, in September 1910 [Hurtig and Kowalle 1988].At the same time, in addition to the pendulum response, other unpredictable, albeit less important, factors can affect the recordings, such as defective decoupling between the two horizon-tal components, or friction between the stylus and the paper; in common practice, unless they evidently affect the recordings, these latter are usually neglected.
The knowledge of the pendulum parameters allows the reconstruction of a theoretical analytical form of the instrument response (Figure 4).However, this is only an approximation of the true response, which can be more or less accurate, and which in general remains unknown and can be very different.In particular, this is the case either at lower frequencies, where small differences in the seismogram amplitude correspond to large differences in the ground displacement, or close to the natural period of the instrument, where the response function is very sensitive to the damping constant h (Figure 5).
Of note, the pendulum damping is usually reported in terms of either the damping ratio f (i.e., the ratio between the amplitude of two succeeding peaks, for free decaying oscillations) or the damping constant h (sometimes also called the damping ratio, representing the ratio of the actual to the critical damping).f and h are related by the equation (Richter [1958], see also Kanamori [1988], for this and further useful information on the instrument response of historical seismometers).Sometimes, there can be confusion as to which parameter is referred to in the bulletins and station books; therefore, some caution should be used in their use.

= -^h
Due to the uncertainties on the actual response of the seismometer described above, in general, the reconstruction of the ground motion by deconvolution of the theoretical response from the original recording is certainly not advisable.Anyway, some tests of the possible effects that uncertainty in the instrument parameters can have on the results of waveform analysis should always be performed.

Source locations
Locations of historical earthquakes were often plagued by poor station distributions or lack of data from existing stations (seismograms and/or bulletins) at the time of the event.Also, arrival-time data can suffer from ON THE RECOVERY AND ANALYSIS OF HISTORICAL SEISMOGRAMS  the presence of large outliers, which might derive from timing errors, erroneous phase identification, or difficulties in reading arrival times on analog, possibly smoked, paper recordings.The presence of outliers, in combination with the use of one-dimensional (1-D) velocitydepth models that can be poor approximations of the real Earth structure, results in inadequately constrained locations, with very large errors.Using regionalized and 3-D models can really help in reducing the uncertainty; however, outliers may still represent a major drawback.
For these reasons, modern standard inversion methods for source location are usually not appropriate for old earthquakes.In recent years, Lomax [2005] developed a location code (NonLinLoc) that is based on the use of an equal-differential-time cost function that gives an efficient and quasi-automatic identification of outliers.The algorithm performs an iterative procedure, which uses an importance-sampling method, which then provides a final location in terms of a probability density function in 3-D space.The Lomax [2005] code has already been used for the study of the sources of historical earthquakes, and it has proven to be very effective in providing well-constrained results.Pino et al. [2008], for instance, applied NonLinLoc to the 1930 Irpinia, (southern Italy) earthquake (M W = 6.6).Here they obtained a reliable hypocenter using a total of 63 readings, with uncertainty comparable to the regional/ teleseismic location of modern events, (38 P and 25 S) from 42 stations.

Earthquake magnitude
Magnitude is one of the fundamental parameters for the characterization of an earthquake source.This was introduced for the first time in 1935 by C. Richter, to determine as quantitatively as possible the strength of an earthquake from the recorded ground motion.Thus, since the first seismographic records, several decades passed before a quantitative estimate of the earthquake strength could be obtained from the information contained in the seismograms.The most common approach to earthquake magnitude determination is to consider it as a function of the logarithm of the maximum value of the A/T ratio, where A is the amplitude of true ground motion and T is the period of the corresponding oscillation; i.e., it depends on the size and period of the recording, and therefore it can be seriously affected by both a poor match between theoretical and actual frequency response functions and possible errors introduced in the waveform processing (e.g., digitization, curvature correction).In particular, a wrong estimate of A/T by a factor of 2 results in a 0.3 error in magnitude estimation.The effect of different values for damping constant h and pendulum free period T 0 on the determination of the magnitude is evident from Figure 5: the largest sensitivity to h is at periods around the natural period of the pendulum, as expected, while magnitudes estimated at wave periods around or larger than the instrument free period T 0 (lower frequencies) can be affected by the assumption of an inappropriate value for T 0 .As is obvious, uncertainties in the magnitude evaluation can be statistically reduced progressively for the increasing numbers of available seismograms [Bath 1981].

Spectral analysis
As stated above, a knowledge of the reliable values for the instrument parameters is crucial for any analysis of seismograms, and sometimes the value reported by original bulletins can be inaccurate.The theoretical response function constructed from the listed constants can be tested, at least qualitatively, by performing spectral analysis of the associated seismogram.
In principle, the spectrum of the pulse of the P-wave Figure 6.Comparison between the spectra of seismograms (green) relative to the recordings of the 1930 Irpinia earthquake at Ebro (Spain) (Mainka horizontal seismographs) and the theoretical spectra computed using different the damping constants h, reported.Small letters indicate the components (north, east).Apparently, the damping values reported in the bulletins -0.40 and 0.28, respectively for the north and the east components -appear to be slightly overestimated.
recorded at teleseismic distances can be considered a good approximation for the spectrum of the seismic moment rate released by the source.Therefore, to assess the reliability of the seismograph theoretical response, the spectrum of the recorded P-wave pulse can be compared to theoretical spectra that are calculated for different source models and are corrected for propagation attenuation and instrument response.This procedure can also be applied to test instrument responses when only regional seismograms are available, and for magnitude (M) <7.0 to 7.5 earthquakes (and source duration < ~20 s), and even for crustal sources.Indeed, at epicentral distances between 8° and 10° and 18° and 20° or beyond, the spectrum of the Pwave train for crustal earthquakes is complicated by the superposition of the PL-wave, upper mantle triplications and depth phases.Typical periods for PL phases are usually larger than 20 s to 30 s, and the depth phases [Lay and Wallace 1995, p. 375], as well as the triplication phases, add amplitude at low frequency, below both the source corner frequency and the instrument natural period for most historical seismographs.
It is worth noting that for event-station couples for which the source duration and instrument natural period are distinct from each other, this kind of spectral analysis can also be a viable technique to check the tabled damping constant (Figure 6), at least qualitatively.In particular, this can be useful when there is no information at all available for h.
Overall, while providing a good check on the instrument parameters, the spectral analysis can also provide information on the characteristics of the earthquake source, such as the seismic moment M 0 and the apparent duration of the source, as demonstrated by Pino et al. [2000] (Figure 7).Alternatively, the seismic moment and, thus, the moment magnitude, can be directly derived from the spectral amplitude, using an analytical source model.For instance, Teves-Costa et al. [1999] and Badal et al. [2000] used this second method (by applying the Brune [1970] source model) to estimate the seismic moment of the 1909 Benavente (Portugal) earthquake, of M W = 6.0, and of some small events (3.64 < < M W < 4.85) that occurred in Spain over the period from 1923 to 1961, from the analysis of historical seismograms recorded at local and regional distances.

Focal mechanism
Today, the source mechanism of modern moderate and large earthquakes is routinely obtained by automatic inversion of three-component, high-quality, digital, broadband waveforms, in real-time or quasireal-time, by matching the data with synthetic Green's ON THE RECOVERY AND ANALYSIS OF HISTORICAL SEISMOGRAMS The seismic moment M 0 obtained by scaling the synthetic spectra to fit the data is also indicated.This kind of test provides a first rough estimate of M 0 (for comparison, Pino et al. [2000] report M 0 = 5.38 × 10 19 Nm as the preferred solution, derived from body waveform inversion); but is also useful as a check for the theoretical instrument parameters.In the synthetic spectra, the peak at 0.067 Hz corresponds to the inverse of the assumed source duration (15 s), while the other peak that is visible in some of the plots, is due to the instrument free period.functions (i.e., moment tensor inversions).Usually, instrument response is deconvolved from the recordings and the horizontal components are rotated along the radial and transversal directions.In these analyses, the waveforms are commonly low-pass filtered (below a few hundredths of a Hz, at least) to reduce the effects of minor lateral structural heterogeneities, thus allowing the adoption of 1-D models for the inner Earth in computing the synthetics.
In principle, these methods can be applied to historical seismograms (see, for instance, Huang et al. [1994]).However, specific limitations arise from the availability and the quality of the data, including the information on the instrument response, as uncertainties in one or more of the instrument responses will seriously affect the performance of the whole inversion.Thus, joint inversion at multiple stations is not always feasible or suitable.For these reasons, for older events, single-station inversions are very frequent [e.g., Huang et al. 1998] and sometimes only attempts with two components and/or a portion of the seismograms can be made.
In general, it would be better not to manipulate the original waveforms, thus applying to the synthetics both the correction for instrument response and the rotation of the horizontal components along the seismograph recording directions, rather than the opposite.Nevertheless, several studies have obtained reliable source mechanisms also by deconvolving the response from digitized data [e.g., Ichinose et al. 2003, Stich et al. 2005].
A simpler method for deriving a source mechanism is the inversion of the first motion polarities, which can be effectively applied to the analysis of historical seismograms [e.g., Pino et al. 2008].In this case, the use of SH-wave polarity, albeit restricted to naturally rotated components, can greatly help in constraining the result (Figure 8).Depending on the relative locations of the available stations and nodal planes, on the focal sphere, sometimes a very few SH-polarity readings can be enough to constrain the solution unequivocally.
Although the inversion of the first motion polarities gives less accurate solutions, sometimes the application of moment tensor inversions to historical seismograms can be problematic, inappropriate or even not possible.Indeed, filtering of historical waveforms at low frequencies can produce very unstable results, which seriously affect the estimated moment tensor.At the same time, the seismograms of the strongest earthquakes, especially when recorded at local or regional distances, can be clipped (i.e., the oscillations induced by the ground motion are greater than the maximum allowed by the seismograph and part of the amplitude is cut out) or truncated (i.e., the recording is stopped because the stylus went off scale and did not returned back) (Figure 9).A further advantage of first motion inversion is that this method can also be applied when the pendulum characteristics are poorly or not known at all.Moreover, as bulletins with indications of the firstarrival polarities are often available, solutions can be attempted for the focal mechanism even without having the original seismograms.
Finally, when the data available are not enough for any inversion, which might be because of number or quality, some information on the seismic source also regarding magnitude and focal depth can be obtained by visually comparing the waveforms with the recordings of recent events with similar locations of the source and receivers, and with known source characteristics [e.g., Kanamori et al. 2010] or through qualitative comparisons between recorded waveforms and synthetics computed for a number of suitable sources with different characteristics [Baroux et al. 2003] (Figure 10).In this second circumstance, if a reliable Earth model is available, the relative amplitude of the different wave trains can be very diagnostic and  Pino et al. [2008] for the 1930 Irpinia (southern Italy) earthquake (M = 6.6), by inversion of the first motion polarity of P-waves and SH body waves.A few SH polarities might be enough to constrain the solution; however, only SH on the naturally rotated (transversal axes of stations located at cardinal directions from the source) components can be used, as artifacts are likely to occur in reconstructing SH motion from the composition of both of the horizontal components, at stations located along other directions.
can help in defining a preferred source mechanism.
As is obvious, all of the above analyses are based on reliable knowledge of the orientation and polarity of the seismograph components, which are usually reported in the station logs or bulletins.No check is possible on the accuracy of these parameters, unless seismograms from the same station are available for an earthquake with a well-known source mechanism.When no reference seismic event can be used, tests can only be performed on the inner consistency among the various components and the seismic ray back-azimuth.

Rupture history
Waveform inversion for deriving the slip distribution with time on the whole fault plane from seismograms recorded at local, regional and teleseismic distances is a well-established technique for the analysis of modern earthquakes.Due to the limitations described above, the application of this method to historical recordings is very problematic.In general, when the data the use of inversion procedures, several constraints must be imposed to have reasonable results, otherwise solutions are likely to be affected by artifacts and/or poor resolution.Several studies have been published describing the applications of these methods to historical seismograms.Batlló et al. [2008] cite and give an extensive description of a large number of these studies, and therefore we will not mention them further here.
More often, the quantity and/or quality of the available recordings do not permit the use of such techniques, and ad-hoc analysis should be designed.In most cases, single stations inversion is the best that can be done, with final comparisons of the results from the useable stations for derivation of reliable source images ON THE RECOVERY AND ANALYSIS OF HISTORICAL SEISMOGRAMS   Baroux et al. [2003]).In some cases, this technique can provide useful and reliable constraints on the source geometry, even when only a single station is available.
(Figure 11).As is obvious, in this case, no inversion can be performed to retrieve a variable slip function on the fault plane, and only an apparent source time function can be obtained.This is usually done from the inversion of body waves, and in particular, on the first-arrival Pwaves.It is worth recalling that for moderate magnitude historical events, in particular, most of the available recordings are at regional distances, due to both the sparse stations locations and the prevalent concentration of observatories in some regions.Under these circumstances, the role of the Earth model used to compute the synthetic seismograms is crucial, and due to the heterogeneity of the upper mantle, it is very difficult to define a structural model that can reproduce the ground motion at regional distances, with the details needed for a reliable determination of the source time function, especially for moderate earthquakes.Thus, the effects of the model on the results should also be checked, to avoid the incorrect interpretation of artifacts that derive from unmodeled structural features (Figure 11).
Despite the basic nature of the technique, wellconstrained and very useful information can still be derived from the results.If the useable stations provide sufficient azimuthal coverage, for instance, the rupture length, orientation, and velocity can be estimated by an investigation of the angular variation of the apparent moment rates (Figure 12), and M 0 can also be evaluated by integrating the source time functions.In the case of earthquakes with only a narrow angle that are spanned by the source-receiver paths (i.e., all of the stations in a restricted area), specific analyses can still be performed to get important indications of the rupture  Pino et al. [2008] inverted ~20 moment rate durations M(t) derived from seismograms of the 1930 Irpinia (southern Italy) earthquake, and derived rupture length, velocity, and direction using a grid search technique, coupled with a simple inversion scheme.In particular, they assumed any possible rupture direction and, for each one, computed the variance of the observed source durations with respect to the expected ones for the best fitting rupture length and velocity (a).The minimum variance corresponds to a rupture occurring for 32 km, toward N100°E, at about 2 km/s.These values provide a good fit to the data (b) and are in good agreement with other independent seismological and geological evidence, supporting the validity of the technique.t), with each one derived from single-station inversion, from Pino et al. [2000] for the 1908 Messina Straits earthquake.For all available stations, these authors also computed two source-time functions (shown with distinct line styles), using two different structural models.Comparisons among all of the results allows the identification of the prominent source features.characteristics.Indeed, considering the equation for the apparent duration t a = (1/o r − cosj r /c i ) [Aki and Richards 2002], exploiting the dependence of t a on both the velocity of the seismic waves considered (c i ; i = P,S) and the angle between the rupture propagation and the source-receiver path j r , the fracture propagation direction can be inferred, and in turn, so can the true moment rate function M(t) (Figure 13).
No matter what method is used to derive M(t), in some cases this function can be used to obtain an approximate picture of the slip distribution along a fault.Indeed, for prevalent unilateral rupture propagation, and constant fracture velocity, the source time function corresponds to the slip integral along the fault width, as a function of the rupture length.Thus, following the method introduced by Kanamori et al. [1992] in their analysis of the 1992 Landers earthquake, a rough approximation of the along-strike slip distribution can be derived by simple algebra (Figure 14).It should be noted that when applied to modern events [e.g., Kanamori et al. 1992, Pino andDi Luccio 2009], for which source models have also been obtained from inversions of high-quality huge datasets, this method has been proven to provide reliable and very useful results, albeit crude and much less resolved than modern standard inversions.
The method described above relies on the determination of moment rate functions.These can be affected by the characteristics of the instruments that are used to convolve the synthetic seismograms, which are different from the actual ones.Indeed, small variations in the pendulum natural period, the damping coefficient, or the paper speed can significantly distort the apparent source time functions, which will produce errors in all of the results that were derived from these moment rate functions (Figure 15).In particular, while no dramatic effects are expected on the source duration, the amplitude of the moment rate that results from the deconvolution of a synthetic Green's function from the recorded waveform is sensitive to small variations in the considered parameters.It should be considered that errors in the seismic moment would map in the M 0 estimate, and in turn, in the amplitude of the retrieved slip function.As shown in Figure 15, deviations as large as those reported by Hurtig and Kowalle [1988] for the Wiechert seismograph in Potsdam (see Section 3) can result in a ~0.2 error in magnitude and a more significant inaccuracy in the estimated slip amplitude.We note that this description is meant to provide only a general picture, as the accuracy of the results obtained from the analysis of historical seismograms should always be checked, case-by-case, using tests on the robustness of the single estimates.

Concluding remarks
For any region, quantitative investigations of earthquakes that occurred in the early stages of the seismological instrument era can provide invaluable information on the modes of stress release in this area.Today, both the theoretical advances and the powerful technological and computational tools made available  Pino et al. [2000] suggested that useful information on the source directivity can be obtained by the comparison of the P-wave and S-wave pulses.In particular, once the P-wave apparent source-time functions were derived, they assumed unilateral propagation and computed the expected S-wave apparent source durations for the two possible rupture propagation directions, along the fault strike, oriented N-S.The synthetics computed for these two hypotheses were then compared with the data, which clearly indicated that the northward direction is highly preferred for the waveform similarities, and also for the consistency of M 0 with those obtained from P-wave modeling (around 5 × 10 19 Nm).  by progressive advances over the last decades have allowed the reappraisal of significant past earthquakes, which can lead to the definition of important source characteristics for these events.
However, due to intrinsic uncertainties, the analysis of historical seismograms presents several added possible sources of error, which derive from both the recording procedures themselves and from the waveform processing.Thus, special care needs to be taken in the handling of these kind of data and also in the interpretation of the results.In the whole procedure, the reliable knowledge that the instrument parameters are as correct as possible has a crucial role.For this reason, as well as recovering the original seismograms, the retrieval of observatory bulletins and station logs is fundamental for allowing any analysis of historical seismograms.
The continuous evolution of seismic instruments and the introduction of long-distance transmission of digital data have meant that the seismic observatories that were distributed in considerable numbers worldwide until a few decades ago are now progressively obsolete.Within this framework, the instruments, data, and documents kept in these structures have rapidly deteriorated, and in many cases they have been destroyed or lost, together with their huge content of information.The invaluable activities being carried out by the projects mentioned above are of great value, as they are specifically aimed at the recovery, preservation, and distribution of the historical seismic data.At the same time, renewed interest in historical seismograms has giving an impulse to the development of specific techniques of processing and analysis.All of these initiatives have great scientific, and more in general, cultural significance, and they deserve increased support, also considering the possible outcomes for seismic-hazard evaluation and risk mitigation.Here, we have tried to give an overview of the whole procedure, from recovering seismograms to the application of modern techniques for retrieving seismic source information, with a focus on the most crucial steps.We have also suggested possible checks for the robustness of the data and for the available instrument characteristics, with a description of the effects of various uncertainties on the results obtained.We believe this provides useful indications for the analysis of historical seismograms, and also for the correct interpretation of the resulting characteristics of the seismic source.

Figure 1 .
Figure 1.(a) The seismogram of the June 7, 1910, Irpinia (southern Italy) earthquake (M = 5.8), recorded at the Cheb (Czech Republic) seismic station using a Belar Zlatorog seismograph.In this optical recording, the signal is too thick compared to the amplitude and period of the first wave-trains, masking the details and making any reading very problematic.(b) The September 7, 1920, Garfagnana (central Italy) seismic event (M = 6.5), as recorded by the Wiechert seismograph operating at Uppsala (Sweden): breaks in the waveform due to the time marks prevent a complete reconstruction of the actual ground motion.(c) The January 13, 1915, Marsica (central Italy) earthquake (M = 7) recorded at the same station.Here, several segments of the waveform are missing, due to the fast oscillation of the arm (from Pino [2011]).

Figure 2 .
Figure 2. Comparison between the original seismogram (left) and the corresponding digitized waveform (right), relative to the 1908 Messina Straits (southern Italy) earthquake (M = 7.1), as recorded by the Wiechert seismograph installed in Plauen (Germany).The original waveform shows distortions produced by the finite length and inclination of the arm.

Figure 3 .
Figure 3.Comparison between two digitized seismograms of the 1909 Lambesc (southern France) earthquake (M = 6.0) derived from the same original recording on a Wiechert instrument in Hamburg (Germany), but processed by different operators.Clear differences in the digitized waveforms are seen, which are even more evident in the low-pass filtered waveforms, showing that data processing might have influences on the final seismograms, and thus on the analysis results.The numbers on the lower plots indicate the filter corner frequency.

Figure 4 .
Figure 4. Examples of response functions for seismographs operating in Europe at the beginning of the 20 th century.Capital and small letters indicate station codes and seismograph component, respectively.The curves represent theoretical reconstruction, based on the instrument characteristics reported in the relevant bulletins.All the instruments were Wiechert seismographs, except for DBN-B (Bosch-Omori) and EBR (Vicentini).

Figure 5 .
Figure5.Response functions obtained by varying the damping constant h (left) and the pendulum free period T 0 (right), starting from initial values of h = 0.456, T 0 = 14 s, and static magnification A = 130.As evident, small uncertainties in these parameters can produce large differences in the estimation of the true ground motion, and in turn, in the evaluation of the earthquake source characteristics.

Figure 7 .
Figure 7.Comparison between spectra of seismograms (red) relative to recordings of the 1908 Messina Straits earthquake and theoretical (green) ones (redrawn after Pino et al. [2000].Station code (capital) and seismograph component (small) are reported on top of each plot.The seismic moment M 0 obtained by scaling the synthetic spectra to fit the data is also indicated.This kind of test provides a first rough estimate of M 0 (for comparison,Pino et al. [2000] report M 0 = 5.38 × 10 19 Nm as the preferred solution, derived from body waveform inversion); but is also useful as a check for the theoretical instrument parameters.In the synthetic spectra, the peak at 0.067 Hz corresponds to the inverse of the assumed source duration (15 s), while the other peak that is visible in some of the plots, is due to the instrument free period.

Figure 8 .
Figure 8.The focal mechanism obtained byPino et al. [2008] for the 1930 Irpinia (southern Italy) earthquake (M = 6.6), by inversion of the first motion polarity of P-waves and SH body waves.A few SH polarities might be enough to constrain the solution; however, only SH on the naturally rotated (transversal axes of stations located at cardinal directions from the source) components can be used, as artifacts are likely to occur in reconstructing SH motion from the composition of both of the horizontal components, at stations located along other directions.

Figure 9 .
Figure 9. (a) Seismogram of the 1930 Irpinia (southern Italy) earthquake, recorded by the Vicentini seismograph installed in Piacenza (northern Italy).The waveform is clipped for by 20 min.(b) Seismogram of the same event at the Firenze (central Italy) Omori-Alfani instrument.The recording is suddenly interrupted because the large oscillation takes the needle off the scale.

Figure 10 .
Figure10.Source geometry estimation for the 1909 Lambesc earthquake, by simple qualitative comparison of data (dashed), as recorded by the Göttingen (Germany) Wiechert vertical and horizontal seismographs, with synthetics (continuous) computed for different focal mechanisms (afterBaroux et al. [2003]).In some cases, this technique can provide useful and reliable constraints on the source geometry, even when only a single station is available.

Figure 12 .
Figure12.Example of determination of fault plane parameters from apparent duration of the source time functions.In this particular case,Pino et al. [2008] inverted ~20 moment rate durations M(t) derived from seismograms of the 1930 Irpinia (southern Italy) earthquake, and derived rupture length, velocity, and direction using a grid search technique, coupled with a simple inversion scheme.In particular, they assumed any possible rupture direction and, for each one, computed the variance of the observed source durations with respect to the expected ones for the best fitting rupture length and velocity (a).The minimum variance corresponds to a rupture occurring for 32 km, toward N100°E, at about 2 km/s.These values provide a good fit to the data (b) and are in good agreement with other independent seismological and geological evidence, supporting the validity of the technique.

Figure 11 .
Figure11.Moment rate functions M(t), with each one derived from single-station inversion, fromPino et al. [2000] for the 1908 Messina Straits earthquake.For all available stations, these authors also computed two source-time functions (shown with distinct line styles), using two different structural models.Comparisons among all of the results allows the identification of the prominent source features.

Figure 13 .
Figure13.Although only a narrow azimuth interval was covered by the available stations for the 1908 Messina Straits earthquake,Pino et al. [2000] suggested that useful information on the source directivity can be obtained by the comparison of the P-wave and S-wave pulses.In particular, once the P-wave apparent source-time functions were derived, they assumed unilateral propagation and computed the expected S-wave apparent source durations for the two possible rupture propagation directions, along the fault strike, oriented N-S.The synthetics computed for these two hypotheses were then compared with the data, which clearly indicated that the northward direction is highly preferred for the waveform similarities, and also for the consistency of M 0 with those obtained from P-wave modeling (around 5 × 10 19 Nm).

Figure 14 .
Figure 14.Along the strike-slip distribution for the 1908 Messina Straits earthquake (continuous line; the shaded area indicates the uncertainty), as derived byPino et al. [2000].These authors applied the simple method introduced byKanamori et al. [1992] for computing the slip distribution from source-time functions.For comparison, slip profiles derived from geodetic modeling performed byBoschi et al. [1989] (dashed-dotted) and De Natale and Pingue[1991]  (dotted) are also shown.The clear similarity between the seismological and geodetic results indicated the effectiveness of the technique.In this specific case, the seismic data allowed the exclusion of the sharp source peak, indicated by one geodetic solution, corresponding to Messina harbor.

Figure 15 .
Figure 15.Inappropriate instrument characteristics affect the resulting source functions.(a) The result of Figure 14 is initially assumed here and several inversions have been performed for varying free period (left), damping (middle), and paper speed (right).The effects on the shape of the moment rate function is quite evident: according to what is shown in Figure 14, these would map in erroneous slip distributions, and as is obvious, (b) on M 0 and M W estimates.Here, circles, diamonds, and triangles correspond respectively to the left, middle, and right panels in (a).