Delineation of Crustal Structure at VLC Seismographic Station Using Joint Inversion of Receiver Func...

Ohaegbuchu H. E., Igboekwe M. U., Chukwu G. U.

Journal of Geosciences and Geomatics OPEN ACCESSPEER-REVIEWED

Delineation of Crustal Structure at VLC Seismographic Station Using Joint Inversion of Receiver Functions and Dispersion Data

Ohaegbuchu H. E.1,, Igboekwe M. U.1, Chukwu G. U.1

1Department of Physics, College of Physical and Applied Sciences, Michael Okpara University of Agriculture, Umudike, Abia State, Nigeria


Joint inversion of body wave receiver functions and dispersion data was used to model the shear wave velocity distribution of the crust and upper mantle below VLC (44.16° N, 10.39° E), a broadband seismographic station in Italy. Receiver functions are primarily sensitive to the shear wave velocity contrasts and vertical travel times, whereas the surface wave dispersion measurements are sensitive to absolute vertical shear wave velocity averages and changes as function of depth. Each data set has inherent lapses but by jointly inverting both we are able to draw on the capabilities of one to compensate the imperfections of the other and this provides better S-wave velocity constraints than we would obtain by inverting either data set individually. The receiver functions were computed from the teleseismic earthquakes recorded by VLC station between 2005 and 2012, while the dispersion curves at regional scale were determined by the Frequency Time Analysis and have been used to obtain tomography maps, using the two-dimensional tomography algorithm developed by Ditmar and Yanovskaya in 1987. The inversion results include a crust with a sharp gradient near the surface (shear velocity changing from 2.15 to 3.4 kms-1 in 5 km) underlain by a 13-km-thick layer with a shear velocity of and 3.4 kms-1 another 15-km- thick layer with a shear velocity of 3.72 kms-1, and an upper mantle with an average shear velocity of 4.4 kms-1. The crust–mantle transition has a significant gradient, with velocity values varying from 3.72 to 4.4 kms-1 at about 32 km depth. This result is also in agreement with shear wave velocity cross-section of the area obtained from ITA-LSO data sets.

Cite this article:

  • Ohaegbuchu H. E., Igboekwe M. U., Chukwu G. U.. Delineation of Crustal Structure at VLC Seismographic Station Using Joint Inversion of Receiver Functions and Dispersion Data. Journal of Geosciences and Geomatics. Vol. 3, No. 4, 2015, pp 109-115.
  • E., Ohaegbuchu H., Igboekwe M. U., and Chukwu G. U.. "Delineation of Crustal Structure at VLC Seismographic Station Using Joint Inversion of Receiver Functions and Dispersion Data." Journal of Geosciences and Geomatics 3.4 (2015): 109-115.
  • E., O. H. , U., I. M. , & U., C. G. (2015). Delineation of Crustal Structure at VLC Seismographic Station Using Joint Inversion of Receiver Functions and Dispersion Data. Journal of Geosciences and Geomatics, 3(4), 109-115.
  • E., Ohaegbuchu H., Igboekwe M. U., and Chukwu G. U.. "Delineation of Crustal Structure at VLC Seismographic Station Using Joint Inversion of Receiver Functions and Dispersion Data." Journal of Geosciences and Geomatics 3, no. 4 (2015): 109-115.

Import into BibTeX Import into EndNote Import into RefMan Import into RefWorks

At a glance: Figures

1. Introduction

Seismic waves are waves of energy that travel through the Earth resulting from earthquakes, explosions or other processes that generate vibration of the Earth. There are two types of seismic waves: body waves and surface waves. Body waves travel through the body of the Earth, they are reflected and transmitted at the interfaces where the seismic velocities change, and obey Snell’s law. Two types of body waves exist, which are compressional waves (P-waves) and shear waves (S-waves). The time delays between the arrival of direct P and converted waves can be used to infer the Earth’s model.

Surface waves travel along the surface of the Earth and are typically the most destructive during an earthquake since they have larger amplitude and longer time duration than body waves. There exist two types of surface waves: Rayleigh wave and Love wave. Surface waves are dispersed in gen-eral which means that different frequencies travel at different velocities. This frequency dependence of velocity carries information about the Earth’s structure and by using numerical techniques, we can obtain S-wave velocity-depth structure of the Earth from a particular dispersion curve.

Receiver functions are time series extracted from three-component seismograms by deconvolving the vertical component waveforms from the radial and tangential components of the seismograms [3, 8]. They show the relative response of the Earth’s structure near the receiver and is a straight forward and simple method of constraining the Earth’s crust and upper mantle structure from teleseismic waveforms recorded at seismic stations.

In 2000, Juli`a et al. proposed a complementary technique for jointly inverting both receiver functions and dispersion data simultaneously to provide a better constraint of the shear wave velocity structure than either data set would individually give. With this technique we have been able to obtain a model of the S-wave velocity distribution with depth beneath VLC station in Italy.

2. Data Preparation and Processing

Here we describe the procedure by which data were prepared and the way we processed them. Joint Inversion method requires two types of input data sets: receiver functions and dispersion curves.

2.1. Receiver Functions: Origin of Data Set

Receiver functions were computed from earthquakes recorded at VLC with coordinates (44.16° N, 10.39° E), a broadband seismographic station in Italy, for a group of signals that sampled the same structure. Figure 1 shows the location of the station as well as other stations within the network. Data set were download from the IRIS website (, and events with magnitude greater than 5.5 and good signal-to-noise ratio were selected (events with SN R ≥ 1.5 were selected while those with SN R ≤ 1.5 were rejected). The time window considered is 2 minutes before P arrival and 10 minutes after P. We chose events within all depth range (deep–shallow) but between 30° and 90° from VLC station. Finally, we selected those events in which the observed and theoretical receiver functions fit at more than 60 percent. Table 1 gives a list of teleseismic earthquakes we used in this study.

Figure 1. Broadband seismographic stations of MEDNET
2.2. Placement of Event Source Parameters into the SAC Header

Having on our hands waveforms recorded by the instrument, then we supply the information necessary for SAC to rotate the horizontal seismograms into the theoretically based radial and tangential directions. Header variables (latitude and longitude, the component azimuth, and the component incident angle) are therefore, set in each of our waveforms, the program doallset.c helps to place the event source parameters into the SAC header and to deconvolve the instrument response from ground velocity in units of m/sec.

2.3. Rotation

Rotation is the second step for preparing receiver functions which will serve as one of the input for the joint inversion. In fact, a station records data in three directions: the vertical Z, North-South N , and East-West E. But they not aligned in the axis of the earthquake and the energy in form of various wave types will be found in each of the recorded components. There are two commonly use rotation system: Z − R − T rotation (2D rotation system) which keep the Z-component still pointing in same direction as in the original ZN E recording and the two horizontal components N and E are rotated into the Radial R and Tangential T components respectively.

Table 1. List of teleseismic records used in the receiver function measurements

L − Q − T (3D rotation system), which rotates the Z component into L component where P-wave energy is concentrated, N component is rotated to the Q component where SV -wave energy is concentrated, and E component is rotated to the T component where SH -wave energy is concentrated. The program doallrot.c helps to do these operations.

2.4. Deconvolution

After rotating the signals, we then computed receiver functions. We first used the program doallp.c to pick and cut the traces for the deconvolution. Then we next used the program doallrf.c to compute receiver functions for three different frequencies 0.5, 1.0 and 2.5Hz. The Gaussian factor α used here is 1.0. Then receiver functions can be visualized using the SAC command ppk, they can also be plotted by ray parameter using plotnps SAC command.

2.5. Surface Wave Dispersion Data, Tomography and Inversion

The surface wave dispersion data were used in this study to investigate the shear velocity structure below VLC station. The following techniques were applied in the sequence: (a) frequency-time analysis, FTAN [1], to measure group velocity dispersion curves of the fundamental mode of Rayleigh waves; (b) two-dimensional tomography [12] to map the distribution of group and phase velocities of Rayleigh waves, plotted on a grid to 10 × 10 compute the cellular dispersion curves; (c) non-linear inversion [6] of the assembled cellular dispersion curves to calculate the set of accepted models for each cell; (d) smoothing optimization algorithms [5] to choose the representative model for each cell and thus to define, for the Italic region, the three-dimensional shear velocity model and its uncertainties.

2.6. The Initial Models

The initial models (Figure 4) of the linear joint inversion procedure are the set of models determined by the non-linear “hedgehog” inversion of the cellular dispersion curves in the Italic region as described by Gonz ́alez et al. [10], which correspond to the cells where the stations are located. In comparison with the LSO solution of hedgehog non-linear inversion (Figure 2) obtained by Brandmayr et al. [4], our initial models (Figure 4) are in agreement since VLC station is in d0 cell. The d0 cell is characterized by a gently thinning crust going from 40 to 30 km thickness. The underlying mantle presents two lithospheric layers, thickening westward down to a depth of about 160 km in cell d0 with VS about 4.70 km/s.

Such models are chosen because they fulfil the following conditions: (1) The stations are within cells, whose models are defined with adequate resolution, located in the region studied by surface wave tomography [10]. (2) These cellular models are the solutions of the non-linear inversion of dispersion curves. (3) In each cell, each model differs from the others by at least ±Pi for one of the free parameters Pi (thickness, VS), where Pi is consistent with the resolving power of the dispersion data, as described by Panza [6]. (4) They allow to minimize the drawbacks intrinsic in the linearization of a non-linear inverse problem.

2.7. The Joint Inversion

For this study, the joint inversion technique was used at VLC seismographic station. The program joint96 accessible in the software package Computer Programs in Seismology [11] inverts simultaneously both surface wave dispersion curves and receiver functions. In the inversion algorithm the damping factor is an important parameter because it balances the tradeoff between resolution and stability [8].

Figure 2. Cellular structural model extended down to 350 km depth for the cell containing VLC station (d0) and its neighbours. Yellow to brown colours represent crustal layers, blue to violet colours indicate mantle layers. Red dots denote all seismic events collected by ISC with magnitude greater than 3 (1904-2006). For each layer VS variability range is reported. The uncertainty on thickness is represented by texture (modified from Brandmayr et al [4])

The damping factor was 0.5 for each iteration. The influence parameter p which controls the relevance or the weight given to the receiver functions or dispersion data as explained in previous chapter, was chosen in such a way that more weight was given to the receiver function rather than to the dispersion data that define the initial models. The joint inversion was limited to the first 40s of the receiver functions because all receiver functions computed have poorly constrained features after 35s (very small amplitude) and our goal is to investigate the crustal and upper mantle layers.

The misfit function controls the variation of the average of the percent of fit between the experimental and the theoretical receiver functions generated by the inversion for all used earthquakes at the considered frequencies.

The iterative process is stopped when the improvement of the misfit function from one iteration to the next is less than 0.05 percent. The inversion procedure does not lead to unique solution due to an intrinsic depth–velocity tradeoff associated with the relative nature of receiver functions [2]. The representative solution was chosen among all models based on the following criteria: (a) the solution has the best percentage of fit for the receiver functions and (b) the solution corresponds to a dispersion curve whose difference with the experimental data at each period is within their corresponding experimental errors and the standard error with respect to observed group velocities [9].

3. Results and Discussions

In this section, we present and discussed main results obtained from the joint inversion of receiver functions and dispersion data in Italy.

3.1. Results

The joint inversion of receiver functions and dispersion data, which we performed for VLC station, allows us to get a relatively fine structure for the model for crust and upper mantle structure in the region. Using the cellular group velocity dispersion curves of VLC station and the teleseismic body wave receiver functions recorded by the station, and jointly inverting both, as displayed in Figure 3, the model of S-wave distribution with depth beneath the station was obtained as shown in Figure 5. Some important features of the lithosphere and asthenosphere system, including the Moho depth, the lithosphere–asthenosphere transition zone and the presence a subducted slab at about 100 km depth have been delineated beneath the seismographic station in Italy. An independent shear wave velocity cross-section (Figure 7) was also obtained from ISC catalogue using an interactive tomographic application of, which shows a significant conformity with our model.

Figure 3. Group velocity dispersion curve (left) and receiver functions (right) jointly inverted. The blue lines are the experimental data while the red lines indicate the chosen theoretical receiver functions, corresponding to the best percent of fit

The depth distribution of seismicity (Figure 8) is used as an additional criterion for appraisal of the cellular models. The revised ISC (2007) catalogue for the period 2005–2012 is used and for the cell which contains VLC station we computed histograms, grouping hypocenters in depth intervals.

The depth grouping considered the uncertainties of the hypocenter’s depth calculation: 5 km for crustal and upper mantle seismicity. The computed histograms are the following: (a) distribution of the logarithm of earthquake number per depth intervals (logN–h); (b) distribution of the logarithm of total energy released by the earthquake per depth intervals (logE–h). Each of these depth distributions is generated for two selections of the events: all events in the considered catalogue and all events with depth not fixed during hypocenter calculations in order to reduce the influence of “standard” depths, as 33 km, in the histogram construction.

Figure 6. Chosen solution for VLC station, result of the joint inversion of receiver functions and dispersion data. The high velocity LID, centered at a depth of about 100 km, seems to reveal the presence of a subducted slab
Figure 7. Shear wave velocity cross-section of the study area with distribution of earthquake hypocenters (black dots) with depth

The logN–h distribution gives the an estimation of the material’s fragility in the relevant depth interval: high earthquake’s frequency is related to more brittleness of materials. The logE–h distribution valuates the strangeness of the material in the relevant depth interval: high energy release is related to high energy accumulation in strong materials. This kind of seismic energy distribution is used by Panza et al. [7] but some strong badly located earthquakes can bias the results.

Figure 8. (a) LogN-h, distribution of number of earthquakes with respect to depth, obtained by grouping hypocenters in 5-km of interval. (b) LogE-h distribution of logarithm of energy released with respect to depth

4. Conclusions

The combined receiver functions and dispersion data provides constraints on the shear velocity of the propagating medium that improves those provided by either data set considered individually, and helps to avoid overinterpretation of each data set. Though this combination may not unambiguously resolve the fine structure of the upper mantle, depending on the bandwidth considered for the dispersion data set. When there is no long period information, independent a priori information must be provided during the inversion process to ensure a reasonable upper mantle in resulting model. When there are no short period dispersion velocities, the upper crust velocity information contained in the main peak of the receiver functions data may not satisfactorily contrain that part of the Earth, and independent contraints may also be needed.


[1]  A. L. Levshin, T. B. Yanovskaya, A. V. Lander, B. G. Bukchin, M. P. Barmin, and L. I. Ratnicova (1989) Surface Seismic Waves in Laterally Inhomogenous Earth. In: Keilis-Borok VI (ed).
In article      PubMed
[2]  C. J. Ammon, G. E. Randall, and G. Zandt (1990) On the nonuniqueness of receiver function inversions. J. Geophys. Res., 95:15303-15318.
In article      View Article
[3]  C. A. Langston (1979) Structure under Mount Rainier, Washington, inferred from the teleseismic body wages. J. Geophys.Res., 84:4749-4762.
In article      View Article
[4]  Enrico Brandmayr, Reneta Blagoeva Raykova, Marco Zuri, Fabio Romanelli, Carlo Doglioni, and Giuliano Francesco Panza (2010) The Lithosphere in Italy: Structure and seismicity. Journal of the Virtual Explorer., 36:20-44.
In article      View Article
[5]  G. Boyandzhiev, E. Brandmayr, T. Pinat, and G. F. Panza (2008) Optimization for nonlinear inversion. Earth and Environmental Science, 19:17-43.
In article      
[6]  G. F. Panza (1981) The resolving power of seismic surface waves with respect to crust and upper mantle structural models. In: Cassinis R (ed) The solution of the inverse problem in geophysical interpretation. Plenum Press New York.
In article      View Article  PubMed
[7]  G. F. Panza, A. Peccerillo, A. Aoudia, and B. M. Farina (2007) Geophysical and petrological modelling of the structure and composition of the crust and upper mantle in complex geodynamic settings: The Tyrrhenian Sea and surroundings. Earth Science Reviews, 80:1-46.
In article      View Article
[8]  J. Juli`a, C. J. Ammon, R. B. Hermann, and A. M. Correig. (2000) Joint inversion of receiver function and surface wave dispersion observation. Geophys. J. Int., 143:99-112.
In article      View Article
[9]  O. Gonz ́alez, L. Alvarez, M. Guidarelli, and G. F. Panza (2007) Crust and upper mantle structure in the Caribbean region by group velocity tomography and regionalization. Pure Appl. Geophys., 164:1985-2007.
In article      
[10]  O’Leary Gonz ́alez, Bladimir Moreno, Fabio Romanelli, and Giuliano F. Panza (2012) Lithospheric structure below seismic stations in Cuba from the joint inversion of Rayleigh surface waves dispersion and receiver functions. Geophys. J. Int., 189:1047-1059.
In article      
[11]  R. B. Hermann and C. J. Ammon (2002) Computer programs in seismology: surface waves, receiver functions and crustal structure.
In article      
[12]  T. B. Yanovskaya and G. M. Molchan (2001) Development of methods for surface wave tomography based on Backus–Gilbert approach. Computational Seismology. In: Keilis-Boro, VI (ed).
In article      
  • CiteULikeCiteULike
  • MendeleyMendeley
  • StumbleUponStumbleUpon
  • Add to DeliciousDelicious
  • FacebookFacebook
  • TwitterTwitter
  • LinkedInLinkedIn