Crustal Thickness of Funafuti Using Receiver Function Analysis

In this study, the characterization of the depth of the Mohorovičić discontinuity under the crust of Funafuti island was determined by analyzing the 3 component seismograms from 54 different earthquake events recorded by the station between 2008 and 2012. These seismograms were from teleseismic earthquakes whose epicenter lay at distances greater than 3000 km from the station. The seismograms were iteratively deconvolved in the time domain to remove the unwanted noise and then stacked to obtain better receiver functions. For analysis of the receiver functions, it was assumed that the range in which the Vp/Vs ratio would lie for the given region would be between 1.60 1.85 and the depth of the discontinuity was assumed to lie between 5 20 km. Analysis of the receiver functions showed that the Mohorovičić discontinuity was at a depth of 11 km and the Vp/Vs ratio was 1.75 for the region.


Introduction
Funafuti is an atoll which is the capital of the island nation Tuvalu.The Tuvalu islands are a group of nine islands, which are located midway between Australia and Hawaii.Historically the nation is believed to be first inhabited by migrated Polynesians.These separate islands were together earlier known as Ellice Islands and were ruled by the British Empire.In 1974 the people of Tuvalu voted for the creation of a separate nation and on 1 st October 1978 their dreams were realized and Tuvalu became fully independent within the Commonwealth.Funafuti was also historically the first site where drilling expeditions were first done and borings down to 400 m were dug by the Royal Society of London to study the coral rocks and reefs.
Receiver function analysis is a very easy and an economical method [1] used to map the subsurface structures.It becomes very expensive if conventional seismic imaging techniques using explosives and geophones are used or if other techniques like wireline logging and downhole logging are used.Receiver function analysis has been used many times in the past [2,3] to find out the depth of the Mohorovičić discontinuity which is the boundary between the lower crust and the upper mantle.

Propagation of Seismic Waves through the Earth
The Earth is stratified into many horizontal layers.Each layer has its characteristic physical and chemical properties.A simple structural model of the Earth is shown in Figure 1.
The Mohorovičić discontinuity (Moho) is defined as the depth at which the velocity of the P waves exceeds 7.6 km/s [4].There is also a sudden increase in the S wave velocity.Typically in continental crust this discontinuity is found at depths between 30 km and 60 km and in oceanic crust it is found between 5 km and 9 km.This sudden change in the P and S wave velocity is attributed to the sudden change in the mineralogical properties [5] of the crustal medium.This sudden change affects the density and elastic properties of the crustal medium which typically control the velocity of the seismic waves.
Receiver functions are composed of teleseismic waves which arrive at a receiver station from very large distance earthquakes.The distances are taken large as the waves can then be approximated to plane wavefronts and the waves then reach the station almost vertically after passing through the Moho [6].At the Moho the P wave splits into the Pp and Ps phases which are recorded by the seismogram.Since the rays are approximated to be near vertical, a simple time difference in arrival times between the phases corresponds to the depth of the discontinuity provided the velocity model is known.The resulting seismogram is a three component seismogram with two horizontal and one vertical components.Theoretically one of the horizontal components must not record any signal, but this does not happen in reality as the Earth is not completely homogenous in structure.The resulting seismograms contain noise in the form source effects, propagation path effects, near receiver function effects and due to the instrument response.There will also be reverberations which will be recorded.The near source effects modify the seismic waves near the source of the signal by creating reverberations and internal reflections due to the geological conditions of the source and they effect the transmission of the teleseismic waves and makes its nature very complex.The propagation path effects are due to the various heterogeneities that the wave encounters when propagating through the Earth which causes it lose energy.The near receiver effects are due to the local geological conditions near the receiver which dominates when the seismic waves approach it.The resultant seismograms are due to the combination of all the above effects as shown Equation ( 1).U V (t), U R (t), U T (t) are the vertical, radial, and the transverse components of the seismogram.The operator * is a convolution operator and t is the time variable.I(t), P(t),

U t S t N P t N I U t S t N P t N I t U t S t N P t N I t
N s , and {N RV , N RR , N RT } are the instrument response effect, path effect, source effect and the corresponding receiver function effects.S(t) is the Earths response function.The resultant seismogram is taken as a convolution of all the effects.To remove all these effects and obtain the Earth's response function at the receiver which is otherwise called as receiver function, a deconvolution operation must be performed.Figure 2 shows a simple illustration of how a receiver function is created after an iterative time domain deconvolution.
Earthquake events used for receiver function analysis of the crustal structure beneath a station must meet several requirements.To ensure that the Pp and Ps wavefields are naturally separated on the vertical and horizontal components of the seismogram, only nearly vertical emergent angled seismic waves must be selected.This is achieved by selecting only those earthquake events which are teleseismic in origin with epicentral distances varying between 30˚ and 90˚ from the receiver.The above approximations are valid when the assumption is made that the travel paths of the up going events progressively becomes vertical as it arrives at the receiver.Further only those earthquakes are chosen whose magnitudes are equal to or greater than 6.0 to ensure a good signal to noise ratio.In addition the seismograms are visually verified to ensure that the P arrivals are impulsive and only then are they sent for further processing.

Deconvolution
Several deconvolution techniques like water-level deconvolution technique [7], iterative time-domain deconvolution technique [8] and the multi spectral correlation techniques have been proposed [9] to calculate the receiver functions.When the signal to noise ratio is high and the signal is wideband, all of the above mentioned techniques are equally good.Thus the best approach to compute receiver functions is to compute it for events of large magnitude.However, for studies which need a good back azimuthal coverage there are never enough observations from all azimuths, and consequently one is forced to incorporate signals from smaller events which leads to instability in the deconvolution operation done by spectral division.Iterative time domain deconvolution technique has been found to work very well even in the presence of noisy data and is also used in receiver function analysis in this study.The principle utilized in the time domain deconvolution technique is that receiver functions when convolved with the vertical component of the seismogram would produce the horizontal component of the seismogram.The of the iterative time domain deconvolution method [10] is least-squares minimization of the difference between the observed horizontal seismogram and a predicted seismogram.This is obtained by the convolution of an iteratively updated spike train and the vertical component of the seismogram.As a first step, the vertical component is cross-correlated with the radial component to estimate the lag of the largest spike in the receiver function.The spike amplitude can be obtained by solving a simple equation as illustrated in [11].The current estimate of the receiver function is convolved with the vertical component of the seismogram and the procedure is repeated to estimate other lags and amplitudes.As new lags and amplitudes appear on the predicted receiver function, the misfit values between the predicted radial component and the observed radial component is reduced.In this method the maximum amplitude spikes are first extracted and then progressively smaller and less significant ones are extracted.The iterations are repeated till an acceptable level of misfit is achieved.The quality of the receiver function is determined by the closeness in the matching between the observed and the computed horizontal component for the corresponding receiver function.The advantage of iterative time domain deconvolution method over spectral division is that it is robust even in the presence of noise and the bandwidth does not seem to have an effect in the deconvolution process.

Stacking of Receiver Functions
The amplitude and the arrival times of the Ps phases generated because of the interfaces vary as a function of back azimuth, the distance between the source and the receiver along the great circle and the event depth.In most of the cases, the seismograms have a considerable amount of noise.One of the main reasons behind stacking is to increase the signal-to-noise ratio in the receiver function.However while stacking it is important to consider the ray parameter.Stacking of receiver functions is performed in narrow bins of 10 degrees of back azimuth [12].Stacking over a narrow bin of back azimuth is bet-ter as it helps in the identification of phases in case of a rapidly varying structure.If R1(t), R2(t),•••Rn(t) are receiver functions, then the stacked receiver function Rs(t) is defined by Equation (2).

Geology of Funafuti
The Funafuti atoll is elongated in a near meridional direction due to it being situated on a volcanic zone [13].
The rock specimens that are commonly found on the atoll [14] are tabulated in Table 1.
The lagoon of Funafuti atoll is the largest in Tuvalu and has a north-south extension of 24.5 km and an east-west extension of 17.5 km.The area of the lagoon is 275 km 2 and extends down to 52 m.The lagoon is submerged with reefs.The cross section of the atoll suggests that there are primarily three foundation rocks [15].The first and the topmost layer is the reef and lagoon calcareous deposits.The second could be firm carbonate rocks or equivalently volcanic rocks.The third layer is mainly oceanic crust and composed with basaltic rocks.The western side of the Funafuti atoll, has larger and bigger pebbles of pumice stones as compared to the eastern side.This is attributed to the Carolinian Monsoon current.There is a West-North-West current that prevails in Funafuti during the trade-wind season and an easterly current during other seasons.The entire coastline of Funafuti is affected by erosion [16].Erosion dominates deposition and the atoll is losing its beaches at a very rapid pace.

Estimating the Crustal Thickness of Funafuti
The principle used in estimating the crustal thickness of Funafuti is outlined in [17].The arrival times of the various phases of the P wave namely Pp, Ps, PpSs, and PsPs was noted.In the naming of the phases except for the first arrival all capitals letters represent down going events and all small letters represent up going events of Table 1.Rock specimens found in Funafuti atoll.

Rock Specimen Properties
Coarse Sand Found in the western sea beach of the atoll and is entirely calcareous.
Calcareous Conglomerate Found at a depth of 10 feet from a bore in the southern sea coast.It is loosely cemented and its composition is similar to beach sand.

Conglomerate
There exists abundant wornout pieces of coral which is very much consolidated with limestone and is heavier than the calcareous conglomerates.
Coral Rock Found abundantly in conglomerates and breccias.
Pumice Pebbles Found in the outer circumferences of the atoll and are well rounded and water worn with a fibrous texture.
the corresponding phases.Except the first two events all the other mentioned events produce a down going amplitude on the seismogram.The PpSs and the PsPs phases are recorded at the same time.Deconvolution of the vertical component with the two horizontal components ensures the creation of the receiver functions from which these phases can be easily observed.The time difference between these observed phases corresponds to the depth of the Moho from which these phases originnated.An initial range of the possible values of the depth and the Vp/Vs ratio for a given region must be assumed.The linear relationship between the observed time differences and these assumed values ensures that a best-fit result of the crustal thickness is obtained.This is achieved by plotting a contour map between Vp/Vs ratio versus the assumed depths.The region where the contours converge corresponds to the depth of the Moho and the Vp/ Vs ratio for the region.The station Funafuti is located at latitude −8.52589˚ and longitude 179.19661˚.More than 54 different earthquake events were analyzed between 2008 and 2012 which occurred between 52.8˚N and 30.3˚N latitude and 132.1˚W and 102.9˚E longitude.The receiver functions were stacked and analysis was done with the assumption that the P wave velocity in the region was 6.0 km/s and the depth of the crust in the region lied between 5 km and 20 km and the Vp/Vs ratio was between 1.60 and 1.80.Figure 3 shows the results of the analysis.
The white region which is also the region where the contours seem to converge corresponds to the depth of the Moho and the Vp/Vs ratio for the given region.The places where the contours seem to touch are due to artifacts and should be ignored.

Conclusion
The depth of the Mohorovičić discontinuity for the Funafuti atoll obtained from receiver function analysis was found to be shallow at 11 km and the Vp/Vs ratio for the given region was around 1.75.

Figure 1 .
Figure 1.A simplified illustration of the internal structural of the Earth.Taken from [1].

Figure 2 .
Figure 2. (a) and (b) show the two horizontal components of the seismogram, and (c) shows the vertical component.Iterative time domain deconvolution gives (d).