Shear Wave Velocity Models Beneath Antarctic Margins Inverted by Genetic Algorithm for Teleseismic Receiver Functions

Study on shear velocity structure of the crust and the uppermost mantle around Antarctic continent was started in the 1960’s by using surface waves of the earthquakes occurred around Antarctic plate (Evison et al., 1960; Kovach and Press, 1961). Permanent seismic stations have been operated from the end of 1980’s at Antarctic margins except for the South Pole (SPA) on the continental ice sheet. A majority of the stations were established as the Federation of Digital Seismographic Networks (FDSN; Butler and Anderson, 2008). The FDSN was composed of several national and governmental organizations such as the Global Seismographic Network (GSN) organized by the Incorporated Research Institutions for Seismology (IRIS), the Australian Government (AG), GEOSCOPE by French, Geo Forschungs Netz (GEOFON) by Germany, MEDNET by Italy, PACIFIC21 (Tsuboi, 1995) by Japan, and others.


Introduction
Study on shear velocity structure of the crust and the uppermost mantle around Antarctic continent was started in the 1960's by using surface waves of the earthquakes occurred around Antarctic plate (Evison et al., 1960;Kovach and Press, 1961).Permanent seismic stations have been operated from the end of 1980's at Antarctic margins except for the South Pole (SPA) on the continental ice sheet.A majority of the stations were established as the Federation of Digital Seismographic Networks (FDSN; Butler and Anderson, 2008).The FDSN was composed of several national and governmental organizations such as the Global Seismographic Network (GSN) organized by the Incorporated Research Institutions for Seismology (IRIS), the Australian Government (AG), GEOSCOPE by French, Geo Forschungs Netz (GEOFON) by Germany, MEDNET by Italy, PACIFIC21 (Tsuboi, 1995) by Japan, and others.
In recent few years, surface wave tomography studies around Antarctic continent and surrounding oceans have been conducted by using the FDSN data (Roult et al., 1994;Ritzwoller et al., 2001;Danesi and Morelli, 2001;Kobayashi and Zhao, 2004).Enderby Land, particularly around the Napier Mountains, was one of the oldest Archaean cratons with a spatial extent about 500 km (Ellis, 1987;Black et al., 1987).However, surface wave analyses could not provide enough spatial resolution for detail discussion about fine crustal structure.Therefore, it is necessary to achieve smaller-scale heterogeneities in the specified area by using recently available broadband waveform data.
In this chapter, precise shear velocity models of the crust and the uppermost mantle were investigated from teleseismic receiver functions beneath several FDSN stations in Antarctica (Fig. 1; MAW; 67.6°S, 62.9°E, SYO; 69.0°S, 39.6°E, DRV; 66.7°S, 140.0°E,VNDA; 77.5°S, 161.9°E,PMSA; 64.8°S, 64.0°W).The obtained velocity models were discussed in relationship with the regional tectonics and crustal evolution of each terrain around the stations.

Methodology and data
The coda parts of teleseismic P-waves contain a significant amount of information on the seismic structure in the vicinity of the recorded stations.The "receiver functions" were defined as the structural response and consisted of P-to-S converted waves and their reverberations, and were most sensitive to the shear velocity beneath the station (Fig. 2).To derive the structural response (receiver functions) beneath the recording station, the sourceequalization method (Langston, 1979) was generally applied for the coda part of teleseismic P-waves.The structural response could be isolated from that of the instrument and effective source function.The followings are the procedure to produce the receiver functions.
Three components (V, R, T) of the total response at a station on a teleseismic P wave are expressed as, where I(t) is the impulse response of the recording system, S(t) is the effective seismic source function, and E(t) implies the impulse response of the earth's structure.
For a steeply incident P wave, source equalization method (Langston, 1979) assumes, ( ) ( ) whereδ(t) is the Dirac delta function.From ( 1) and ( 2), Thus the observed structual response E R (t) and E T (t) are available by deconvolving D V (t) from D R (t) and D T (t).Since the receiver functions are sensitive to P-to-S conversions through the interfaces beneath the recording station, the waveform inversion result could produce a shear velocity structure (Owens et al., 1984;Kind et al., 1995).By applying these conventional procedures, receiver functions were obtained at five permanent FDSN stations in Antarctica.

Ray diagram Receiver function
Here, we introduce the actual procedure to create receiver functions as an example for the station MAW.Before the inversion, we used the stacked receiver functions for all 20 radial components in the backazimuth group within 70°.The incoherent noise could be suppressed by stacking the waveforms, while the coherent signals were enhanced.Dataset of teleseismic waveforms for the other four stations were the same as used in the linearlized inversion analyses by Kanao et al. (2002).
Inversion of the receiver functions to recover crustal and uppermost mantle structure have been widely recognized as sensitive to the starting model if a conventional linearization scheme was employed (Ammon et al., 1990).Prior to this study, a linearlized time domain inversion was applied to determine the velocity model for several Antarctic stations (Kanao, 1997;Kanao et al., 2002).The starting model dependency might be excluded by employing a non-linear inversion scheme based on a Genetic Algorithm (GA; Sambridge and Drijkoningen, 1992;Shibutani et al., 1996).

What is genetic algorithm ?
GA in non-linear optimization include three steps: The GA approach makes use of a 'cloud' or 'population' of the models to minimize the dependence on a starting model; a set of 'biological' analogues are used to produce new generations of the models from previous generations, with preferential development of the models with a good fit between observed and theoretical receiver functions.Figure 3  In this study, a non-linear GA was applied for the stacked radial receiver functions of each station.Figure 4 represented a flow chart of GA for receiver function inversion.Beginning with a randomly generated initial population and corresponding misfit values, which were defined by the square sum of the difference between the receiver function predicted for each model and that obtained from observed waveforms, succeeding populations were created by selection, crossover and mutation procedures.
The approach provided a good sampling of the model space, and enabled the estimation of the shear-wave speed distribution within the crust, along with an indication of the ratio between Vp and Vs.Many models with an acceptable fit to the data were generated during the inversion, and a stable crustal model was produced by employing a weighted average of the best 1,000 models encountered in the development of GA.The weighting criterion was based on the inverse of the misfit for each model, so that the best-fitting model could have the greatest influence.
For the inversion, the crust and the uppermost mantle down to 60 km were assumed to be composed from five major layers (Table 1).The model parameters in each layer were the thickness, Vs at the upper and the lower boundaries, in addition to the Vp/Vs ratio.The Vs for each layer was constructed by linearly connecting the values at the upper and the lower boundaries, so as to give a sequence of constant velocity-gradient segments separated by velocity discontinuities.The thickness and the upper and the lower limit in each layer were defined after the averaged continental crust.'Q α 'and 'Q β ' values were assumed to be fixed in each layer on the basis of Coda-Q inversion results after Kanao and Akamatsu (1995).A smoothness constraint in the inversion was implemented by minimizing a roughness norm of the velocity model (Ammon et al., 1990).The thickness and the upper and the lower limit in each layer were defined after the averaged continental crust.'Q α 'and 'Q β ' were assumed to be fixed in each layer by referring from Coda-Q inversion results after Kanao and Akamatsu (1995) After examining the trade-off curves between the model roughness and waveform-fit residuals, we selected the most suitable pair of the above parameters.A number of iterations up to 200 times were conducted in the inversion in order to reduce the waveform-fit residuals (misfit-values) to an acceptable value, and the most stable solutions were adopted as the final models (Fig. 5).We obtained 50 population models for the every iteration.In total, we selected 10,000 models to determine the best fitted.

Basement
The waveform fits between synthetic and observed receiver functions were generally adopted, implying the adequate inversion procedures with reasonable smoothness constrained.Figure 6 represented the synthetic radial receiver functions by assuming the S-wave models and the Vp/Vs ratio determined by the averaged one for the best 1,000 models in GA inversion (broken traces), compared with the observed mean (upper solid trace) and +/-1 standard error bounding (lower two solid traces) of the stacked receiver functions.
There were several noticeable later phases for all traces after the P-arrival.For example, large amplitudes were identified around 4-5 s, which were considered to be the directly converted Ps at the Moho discontinuity.Intra-crustal converted phases were identified around 1-2 s and 2.5-3.5 s, implying the mid-crustal velocity discontinuities.Later phases, after around 7 s, had a rather worse waveform fitting compared with the earlier phases, because of relatively poor signal-to-noise ratios for these later phases.

Results and discussion
In this section, we discussed the resultant shear velocity models for the individual FDSN station.Here, the averaged Vs models for the best 1000 misfits in GA inversion were mainly discussed (red lines in Figs.7a, 7b and 7c).The inverted velocity model beneath MAW had a very sharp Moho discontinuity at a depth of 44 km (Fig. 7a).A discontinuity between the middle and the lower crust were recognized at 34 km depth.In Mac.Robertson Land, where MAW is located (Fig. 1), late-Proterozoic metamorphic events generated the granulite facies rocks in upper part of the crust (Tingey, 1982;Sheraton et al., 1987).The sharp and fairly deeper Moho around 44 km depth might have a relationship with the metamorphism of the Rayner Complex besides the Archaean Napier Complex.The intrusion of charnockites around MAW was an evidence of the compression tectonic setting in the Proterozoic mobile belt (Young and Ellis, 1991).Depletions of heavy rare-earth elements in the low-Ti charnockites suggested that garnet was a residual phase in partial melting, which required high pressures and an overthickened crust.The deeper crustal thickness obtained from GA inversion at MAW appeared to be correspond with a signature of the crustal root what now have been remained as the deepened architecture comparing with the adjacent areas in Enderby Land.

MAW S-velocity model
The resultant velocity model around DRV (Fig. 7b, left) indicated a fairly sharp Moho at depths of 28 km.A high-velocity zone appeared in the upper and the lower crustal depths.
A relatively lower velocity zone was obtained at depths around 20 km, which lied between the above two high-velocity zones.The velocities of the topmost mantle had lower values of about 4.2 km/s.In Adelie Land, where DRV is belonging (Fig. 1), a metamorphic event occurred in early-Proterozoic age (Bellair and Delbos, 1962).A rather sharp Moho and fluctuations of the crustal velocities might had been developed during the metamorphic event of the Adelie Land.In addition, high velocity zones in the upper crust together with a low-velocity discontinuity in the middle crust might be related to the early-Proterozoic tectonothermal activities.
A sharp Moho discontinuity was determined approximately at 40 km beneath SYO (Fig. 7b, right), in the Lüzow-Holm Bay region.The Moho depth was consistent with that obtained from previous large scale deep refraction / wide-angle reflection surveys around the region (Ikami and Ito, 1986;Yoshii et al., 2004).Velocity jumps were identified at 12 km and 20 km depths, which corresponded to the discontinuities between the upper-middle crust and middle-lower crust, respectively.The latter discontinuity between the middle and the lower crust significantly coincided with the depths by the recent refraction / wide-angle reflection study around the SYO (Miyamachi et al., 2003).High velocity zones in the upper crust were presumably corresponding to the granulite facies metamorphic rocks appeared in surface geology.The obtained velocities in the upper part of the crust were consistent with the velocities of granulite facies rocks found from high-pressure laboratory measurements (Christensen and Mooney, 1995;Kitamura et al., 2001).The considerable crustal evolution models to explain the velocity variations within the crust might be related to the compressional stress during the early-Paleozoic metamorphism in the Lützow-Holm Bay region (Hiroi et al., 1991;Shiraishi et al., 1994).
As for the Antarctic Peninsular, very broad Moho discontinuity was found around 36 km depths beneath PMSA (Fig. 7c, left).Several previously conducted refraction / wide-angle reflection experiments had revealed a complicated Moho topography around the region (Sroda et al., 1997;Grad et al., 2002).They determined the thickness of the crust in the range of 36-42 km at coastal area of the Peninsula, in contrast, decreased to about 25-28 km toward the Pacific Ocean.The dipping Moho obtained from our results supported a possibility of the transition zone between the oceanic and continental crust in the Antarctic Peninsula.Since the GA inversion applied for this study had assumed a uniformed structure composed of the flat-lying layers, the dipping structure cannot be directly identified.However, the obtained crust-mantle boundary was considered to reflect the averaged structure for the complicated regime in the vicinity of the Moho discontinuity.

S-velocity model : DRV, SYO
Broadening low velocity zones around 30 km depths and transitional Moho with few km widths at VNDA (Fig. 7c, right) might be involved in the uplift mechanism on the West Antarctic Rift System (WARS) nearby the Trans Antarctic Mountains (TAM) (Smith and Drewry, 1984;Stern and ten Brink, 1989;Ten Brink et al., 1997).Around station VNDA in the Terra Nova Bay region, the Moho depth was already estimated from Ps converted phases of the receiver functions by temporary seismic array data (Di Bona et al., 1997).They pointed out a thinned crust with thickness drastically varied from 17 to 29 km, which implied a transitional zone between East and West Antarctica, which crossing the WARS.The other seismic refraction data indicated the same regime of the Moho depths involving the crustal rift system at TAM (Vedova et al., 1997).Wiens et al. (2003;2006) conducted the TransAntarctic Mountains SEISmic experiment (TAMSEIS) around the region and obtained a detailed distribution of the crustal thickness by receiver function analyses (Lawrence et al., 2006).Their results also indicated relatively shallow Moho depths together with low velocity zones beneath VNDA.
Figure 8 demonstrated a comparison of the Moho depths by three different methods of seismic refraction studies, gravity-based estimates (after Von Frese et al., 1999), together with receiver function GA inversion by this study.In spite of the existence of small differences in estimating the Moho depths between three methods, it might be mentioned that a principal difference between the East and West Antarctica, as well as the Antarctic Peninsula, was remarkably identified.In order to establish a crustal model of the whole regions in Antarctica, available broadband waveform data of the other seismic stations, such as SPA (90.0°S),CSY (66.3°S, 110.5°E),SBA (77.8°S, 166.8°E) and the other temporary stations should be compiled for comparison.
During the International Polar Year (IPY 2007(IPY -2008)), a major geo-science program had been conducted such as the 'Antarctica`s GAmburtsev Province / GAmburtsev Mountain SEISmic experiment (AGAP/GAMSEIS)' (Wiens, 2007).The AGAP/GAMSEIS was an internationally coordinated deployment with few tens of broadband seismographs over the huge area of continental ice sheet in East Antarctica.The investigations on the high plateau inside the ice covered continent could surely provide detail information on the crustal thickness and mantle structure (Hansen et al., 2010).In contrast, the 'Polar Earth Observing Network (POLENET; http://www.polenet.org/)'was another major contribution to the IPY by establishing a geophysical network mostly weighted on West Antarctica.
The accumulated seismic data during the IPY will be utilized to clarify the heterogeneous structure of the crust and upper mantle, as well as the Earth's deep interior, including the features such as the Core-Mantle-Boundary (CMB) and the lowermost mantle layer (D" zone).The broadband seismic arrays in the Antarctic at IPY and beyond have been conducting a significant contribution to FDSN as viewed from high southern latitude.Mapping of the crustal velocities beneath the whole Antarctic continent could firmly address for the advance in interpreting the difference of structure and tectonics in various terrains of Gondwana super-continent.

Conclusions
In this chapter, seismic shear velocity models of the crust and the uppermost mantle were investigated by teleseismic receiver functions beneath the FDSN stations in Antarctica.In order to eliminate the starting model dependency, a non-linear GA was adopted in time domain inversion of the receiver functions.The shear velocity model beneath MAW represented a sharp Moho boundary at 44 km depth.A fairly sharp Moho was identified around 28 km and 40 km depths beneath DRV and SYO, respectively.Shear velocity variations in the crust for these stations might have a relationship with the lithologic variations of metamorphic rocks in the shallow crustal depths.Broadening low-velocity zones around 30 km depths with transitional crust-mantle boundary were identified at VNDA; which might be involved in the WARS associated with elevation of TAM.Moreover, a broad crust-mantle transition was determined around PMSA, in Antarctic Peninsula.These variations in shear velocities within the crust presumably reflected the tectonic history of each terrain where these permanent stations are belonging.

Acknowledgement
The authors express sincerely thankfulness to the data management centers for all the FDSN related stations in Antarctica and collaborators who were involved in data acquisition of the broadband seismic waveforms.We would like to acknowledge IRIS/Data Management System (DMS) for their grateful assistance to get PMSA and VNDA data accessible.Australian Government (AG) and the Australian Seismological Center of the Australian Geological Survey Organization (AGSO) took a grateful effort to make MAW data accessible.They also wish to acknowledge to the Ecole et Observatoire de Physique du Globe for their useful information on GEOSCOPE's DRV data accessibility.Data of SYO were available from Polar Data Center of NIPR.

Fig. 1 .
Fig. 1.Map showing the location of analyzed seismic stations by GA inversion (Solid triangles; SYO, MAW, DRV, PMSA and VNDA) in Antarctic continental margin and the related regional localities.Open triangle stations (CSY, SPA and SBA) are planned to be analyzed in future.Solid red circles are stations raveled in Von Frese et al. (1999).Alphabet numerals are location in Antarctica raveled after Von Frese et al., (1999), the same representation as in Fig. 8

Fig. 2 .
Fig. 2. (upper) The receiver functions (RF, the crustal response) are consist of P-to-S converted waves and their reverberations, which are most sensitive to the shear velocity beneath the station.(lower) the observed receiver functions (E R (t) and E T (t)) can be obtained by deconvolving the original waveforms D V (t) from D R (t) and D T (t)

FlowFig. 4 .
Fig. 4. Flow chart of Genetic Algorithm (GA) for receiver function inversion.Beginning with a randomly generated initial population and corresponding misfit values which are defined by square sum of the difference between the receiver function predicted for each model and that obtained from observed waveforms, succeeding populations are created by selection, crossover and mutation

Fig. 5 .
Fig. 5. Misfit values vs. the number of iteration during the GA inversion for an example of MAW.Variations in the mean, the minimum and the maximum misfit values for each population are drawn to be reached into the stable values Fig. 7a.Seismic S-velocity models for MAW by GA inversion.For the S-wave velocity, all 10,000 models searched in GA inversion are shown as the light gray shaded area.The best 1,000 models are represented by the yellow to green area.The best model and the averaged model for the best 1,000 are shown by the blue line and the red line, respectively.For the Vp/Vs ratio, the light blue solid line corresponds to the averaged model Fig. 7b.Seismic S velocity models for DRV and SYO by GA inversion

Fig. 8 .
Fig. 8.Comparison of Moho depths determined from seismic refraction studies (open dots with points in center) and gravity-based estimates (open squares) (Von Frese et al., 1999), together with receiver function inversion study for permanent seismic stations (Solid diamonds; SYO, MAW, DRV, PMSA and VNDA) by this study.Alphabet numerals are location in Antarctica raveled after Von Frese et al., (1999), the same representation as in Fig. 1

Table 1 .
Model parameters in GA receiver function inversion.'Vs (upper)' and 'Vs (lower)' are the S wave velocity at the upper and the lower boundaries of each layer.The 'Lower' and the 'upper' for four variables indicate the lower and the upper bounds.
Shear Wave Velocity Models Beneath Antarctic Margins Inverted by Genetic Algorithm for Teleseismic Receiver Functions 249 www.intechopen.com