The velocity structures derived from anisotropic inversion of teleseismic receiver functions around Flores Island, Indonesia

We determine the shear wave velocity structure beneath a seismic station located in the Flores Island using anisotropic inversion of Neighbourhood algorithm. The study area is located in the boundary transition zone from subduction to collision in eastern part of Indonesia, providing an ideal and unique place for crustal anisotropic studies. The previous crustal study using receiver function analysis conducted around this region found the amplitude and timing variations of arrivals as a function of back azimuth. These variations suggest the presence of anisotropic structures underneath the region. Thus, in this study we reanalyze the velocity structure derived from receiver functions by incorporating anisotropic parameters in the inversion process. In general, the best fit model derived from the anisotropic inversion is comparable with the previous study. The crustal thickness beneath the study area is around 38 km and the model also shows P and S waves anisotropic layer at the upper and lower crust. The difference of fast anisotropy axis observed on these layers may reveal the different mechanism of deformation between the shallow and deeper parts of the crust underneath the study area.


Introduction
A receiver function is a time series waveform derived from three-component seismograms providing information about the Earth structure underneath a seismic station. The waveform contains P to S converted waves due to the presence of seismic discontinuity underneath the seismic station. In receiver function modelling, if the earth structure is flat, homogenous and isotropic, the receiver function in transverse component provides zero energy. To provide the Earth structure information, receiver function is usually analyzed by inverting the waveform based on assumption of 1-D isotropic earth structure [i.e. 1,2]. This technique is appropriate for simple earth structures involving small numbers of receiver function records. However, the real Earth is very complex containing heterogeneous structures (i.e. dipping interfaces or anisotropic body), thus the tangential component of receiver function will exhibit significant energy [3]. The heterogeneous and anisotropic behaviour of the earth may also be identified by the radial component response showing variation in amplitude of the PS phase as a function of back azimuth. Therefore, involving the complex-structure parameters on the inversion of receiver functions is very important in order to provide better image of the earth structure.
The study area is situated at the boundary between Sunda arc to the west and Banda arc to the east (figure 1). This arc boundary is very complex since it experiences tectonic transition from subduction of Indo-Australian oceanic plate to collision of Indo-Australian continental plate. The convergence rates vary from about 70 mm/yr in the subduction part to 15 mm/yr along the collision part [4,5]. The  [6,7]. Previous geophysical study using shear wave splitting analysis observed crustal anisotropic behaviour related to the dynamic process of the study area [8]. They observed that the anisotropic behaviour in this region is not only governed by stress field, but also caused by the contribution of complex structural fabrics. Furthermore, the receiver function study previously conducted in this area also found anisotropic pattern on the receiver function waveforms [2]. They observed that the receiver functions display a clear amplitude variation of PS converted phases with back azimuth indicating the presence of anisotropic body or dipping interface underneath the seismic station. Although those studies can reveal general anisotropic behaviour and crustal thickness underneath the station, the information about the detail characteristic of complex crustal structure around the study area is remain unclear. Therefore, in this study we provide new perspective about crustal characterization by using complex structure parameters on the inversion of receiver functions.

Method
In this work, we determine the shear wave velocity structure around Flores Island using the same receiver function dataset computed by the previous work [2]. The receiver functions are extracted from seismic waveforms of teleseismic events with magnitude greater than 6 (figure 1). The seismic waveforms are inspected manually to choose the best quality waveform with clear P-arrivals. The waveforms are corrected from the mean, linear trend and effect of instrument response before computing receiver functions, the seismograms are then rotated into radial and transverse components. Finally, the receiver functions are obtained by deconvolving the vertical component of teleseismic waveform from both corresponding radial and transverse components using the iterative deconvolution method [9]. The stable receiver functions in this study are obtained by applying a Gaussian filter with a width factor of 1.5. This width factor effectively removes the frequency noise above 0.75 Hz. Only high quality receiver functions are then used for further analysis by performing at least a 90 % misfit of the radial component. In total, we have 71 high quality receiver functions used in this analysis (figure 2).  The next step for this analysis is inversion process to obtain the velocity structure with complex structure parameters. Since this step requires all traces for the computation process, the computing time will become a big problem. Thus, to resolve this issue, the inversion process is then conducted by two inversion processes. First, we apply 1-D isotropic case to constraint the earth structure pattern and seismic velocity. The second inversion, we then perform 3-D complex structure inversion by considering dipping and anisotropic parameters with the initial model of 1-D isotropic velocity structure derived from the previous inversion. Both velocity structures from the inversion processes are produced using the non-linear neighbourhood inversion method developed by Frederiksen et al. in [10]. We run a number of trials to set the appropriate tuneable parameters with 1000 iterations and a sample size of 20 models. Each inversion process then produces 20020 velocity models. The best inversion solution is that which provides the least misfit between the observed and calculated receiver functions.

Result and discussion
For isotropic inversion, we consider receiver functions with clear PS arrival and back azimuth ranging from 84° to 134°. This back azimuth selection is chosen as the procedure assumes that 1-D isotropic profile should not generate back azimuth variation of receiver functions amplitude at PS arrival. We also only use receiver functions at radial component, as receiver fucntions at the transverse component are assumed to have zero energy for the isotropic case. Our 1-D inversion result shows a good fit between observed and predicted data revealing some major features (figure 3). The velocity profile obtained from this inversion is also comparable to that obtained from the previous study [2] as shown in figure 4. For example, both models show Moho layer at almost similar depth.  For the second step, we add anisotropic and interface dip parameters in the inversion process using the initial model obtained from the first inversion. The inversion requires modelling for all receiver function data set on the two components radial and transverse causing very time consuming on the computational process. In general, the inversion result with minimal misfit shows good fitting for both components ( figure 3 and figure 4) first and Moho arrivals, and also shows quite clear amplitude variation of PS phase compared to that shown by the observed model. This discrepancy may be due to the effect of phase arrivals from dipping structure below the lower crust such as subducting slab. These arrivals are found as negative signals followed by strong positive signals at 7-10 s on the observed radial component. Since we do not incorporate parameters of this dipping structure on the inversion process, thus we limit our analysis for anisotropic and dipping structure located in the crust. For the transverse component, the predicted model recovers quit well for the back azimuthal variation of the significant energy at 0-1 s. Figure 4. the minimal misfit of number of models derived from the inversion process. Figure 5. The best fitting S velocity profile resulted from isotropic and 3-D anisotropic inversion, compared also with the previous model obtained by Syuhada et al. in [2]. The strike, dip and anisotropy characteristics of the best model is also shown in this profile. The best fitting velocity model resulted from the second inversion involves thinner sediment layer compared to that obtained from the first inversion (figure 5). In the next layer, we identify the presence of dipping anisotropic structure with 7 km thick. This structure has P and S anisotropy around 4.6% and 2.6%, respectively, with the fast axis trending 57°; while the interface strikes at S52°E and dips 14°. Another dipping anisotropic layer also found in the lower crust with P and S anisotropy about 7.9% and 3%, respectively, with northward fast direction plunging at 31°. This layer has 12.3 km thick striking S60°E and dipping -7°. Our result is quite interesting as the fast direction of anisotropic layer in the upper crust is almost perpendicular to the fast orientation found in the lower crust. The result implies that the deformation process in the shallow and lower part of the crust beneath the region may have different mechanism. In the upper crust, the cause anisotropy is usually controlled by stress induced anisotropy as also suggested, for example, by Crampin in [11]. Under the influence of tectonic stress, the fast direction of anisotropic layer will be parallel to the stress and strain direction. The fast direction of anisotropic body found in the upper crust underneath the seismic station is generally parallel with the axis of principal compressional strain rate computed by Bock et al. in [5]. In other hand, the fast orientation of anisotropic layer found in the lower crust of the region is perpendicular to the maximum horizontal compressional axis of the strain rate. This can be explained that Flores Island may be under crustal thickening process as consequence of the insertion of Australian continental crust underneath the Banda arc as proposed by other studies [12,13,14]. This interaction causes friction in the plate boundary producing a weak or fracture zone in the lower crust. The presence of fluid as suggested by previous receiver function study then keeps the fracture to be open under high pressure in the lower crust and provides the fast orientation to be perpendicular to the principal stress [8].
In addition, the percent anisotropic values obtained from this inversion are roughly comparable with theoretical crustal anisotropy calculated, for example, by Chang et al. in [15]. Crampin in [11] and Crampin and Booth in [16] also proposed that the typical percentage of crustal anisotropy ranges from 0.5% to 10%. Another interesting feature from this profile is that the crustal-mantle boundary layer is observed at around 38 km depth. This estimation is comparable with the Moho depth estimation obtained from the previous receiver function study [2]. This study suggested that the crustal thickness underneath Flores Island is mainly governed by the interaction of the Australian continental part inserting underneath southern part of Banda arc. This insertion may increase the crustal thickness around Flores Island [13,14].

Conclusions
In this work, we successfully model dipping crustal anisotropic structure beneath a seismic station located in the Flores Island using receiver function analysis. In this model, at least two dipping anisotropic layers have been identified located in the upper and lower crust of the study area. The model also shows that the both layers have different fast direction of anisotropy. This result suggests that the both layers may be deformed with different mechanism. This behaviour is possibly as result of the interaction of Australian continental crust inserting underneath the southern part of Banda arc. In general, the model obtained from this study is also in good agreement with the geophysical studies previously conducted around the region.