Shear wave splitting of the 2018 Lombok earthquake aftershock area, Indonesia

Lombok is one of the islands in the transitional zone from the Sunda Arc to the Banda Arc, Indonesia. In the mid-2018, the island of Lombok was shaken by a series of strong earthquakes, started with a moment magnitude (Mw) 6.4 earthquake on July 29, 2018 followed by earthquakes on August 5 (Mw 7.0), August 9 (Mw 5.9), and August 19 (Mw 6.3 and 6.9). Some researchers suggested that this phenomenon occurred due to a segmentation rupture in the northern part of Lombok Island. This study aims to obtain information on the distribution of the Lombok earthquake fault zone 2018 and also to understand the character of seismic anisotropy around the Lombok earthquake fault zone 2018 through Shear Wave Splitting (SWS) study. Splitting, or S-wave separation, occurs when the S wave passes through an anisotropic medium. The S wave is split into fast and slow S waves with almost orthogonal polarizations and has parameters such as delay time and polarization direction of the fast S wave. To determine the SWS parameters, we used a Lombok earthquake aftershock data set recorded from 4 August to 9 September 2018, using 16 seismographic stations. The steps taken to obtain the SWS parameters are event selection, windowing using short time Fourier transform, and rotation-correlation process. The results of the SWS analysis indicate that the fast polarization directions probably have a linkage with the local fault system and the fault related to the Lombok earthquake fault zone.


Introduction
In an anisotropic medium, shear waves tend to be polarized into approximately orthogonal fast (S fast ) and slow (S slow ) directions (Crampin et al. 1978). The phenomenon of shear wave splitting (SWS) can be characterized by measuring the time delay between the two polarized shear waves (δt) and S fast polarization direction (ф). In SWS analysis, ф represents the polarization direction of the first arriving S wave, which can be viewed by the anisotropic symmetry axes orientation and the propagation direction (Savage et al. 2016). Meanwhile, δt is indicated by the time difference between the fast and slow shear wave waveforms, and is a non-linear product that depends on the anisotropic strength and path length through the anisotropic medium (Savage et al. 2016). SWS phenomenon is generated by anisotropic mediums. The anisotropy condition in a medium is defined by Sasmi et al. Geoscience Letters (2023) 10:7 the dependence of seismic velocity on a certain direction when the seismic wave travels through an anisotropic medium. There are two forms of anisotropy: extrinsic and intrinsic anisotropy. Extrinsic anisotropy is well known for Shape Preferred Orientation (SPO) as it is related to the alignment of the structural features inside the earth and is mainly observed in the crust. Extrinsic anisotropy form is created either due to stress or structure-induced anisotropy. For the second form of seismic anisotropy, intrinsic anisotropy is mostly produced by lattice preferred orientation (LPO) or intrinsic orientation of the anisotropic minerals (Kendall 2000). This anisotropy form is the dominant anisotropic source in the upper mantle. In the crust, anisotropic mineral alignment forming foliation planes produces an anisotropic source through structure-induced anisotropy.
As discussed previously (Kendall 2000), anisotropy in the crust is dominated by extrinsic anisotropy through stress-induced and structure-induced anisotropy. Stressinduced anisotropy can be described as the preferential orientation of the open microcracks under the influence of an active stress field (Boness and Zoback 2006). Thus, the pattern of compressional strain direction in a region may affect the preferred direction of local microcracks. Consequently, stress-induced anisotropy controls the S fast to propagate in a similar direction to the compressional strain of the area. At deeper crust, the effect of this anisotropy usually decreases with increasing the confining pressure closing the microcrack in all directions. However, the presence of fluid-filled microcracks may prevent the cracks from closing contributing to crustal anisotropy (e.g., Hiramatsu and Iidaka 2015). For the second mechanism, structure-induced anisotropy is seismic anisotropy related to aligned macroscopic fabrics such as fault or strike slip filled with fluid (Zhang et al. 2007;Tadokoro and Ando 1999), foliation plane of anisotropic rocks, and sedimentary bedding system (Valcke et al. 2006;Kern and Wenk 1990). In this case, structure-induced anisotropy causes S fast to be polarized parallel to the main strike or foliation plane of the anisotropic medium. SWS analysis has been applied for seismic anisotropy studies either in the crust or mantle (e.g., Audet 2014; Boness and Zoback 2006). In the crustal regime, SWS is usually generated by the presence of local faults or sedimentary layering systems, and observable from near field seismic stations (e.g., Chang et al. 2009). Meanwhile, alignment of anisotropic minerals deformed through dislocation creep may produce SWS phenomenon in the mantle region. In this research, SWS analysis was applied for studying crustal anisotropy in Lombok Island region, Indonesia, using local earthquake waveforms extracted from the local temporary seismic stations.
Lombok Island is located near the border area between the Sunda Arc and the Banda Arc tectonic systems, Indonesia (Fig. 1). The seismicity around Lombok Island is mainly governed by two tectonic arrangements around the area: the subduction zone in the south (Hamilton 1977;1979), and the back arc thrust system in the north (Hamilton 1977(Hamilton , 1979Silver et al. 1983). Some researchers also suggest a subduction structure in the north of Lombok Island called the Flores Oceanic Crust/FOC (McCaffrey and Nabelek 1984;Widiyantoro et al. 2011;Zubaidah et al. 2014;Sasmi et al. 2020;Afif et al. 2021). In addition, Lombok Island is flanked by the Lombok Strait Strike Slip Fault in the west, and the Sumbawa Strait Shear Fault (Sumbawa Strait Strike Slip Fault) in the east (Irsyam et al. 2017). On July 28 2018 22:47:37 UʈTC an earthquake of moment magnitude (M w ) 6.4 hit Lombok Island (Badan Meteorologi, Klimatologi Geofisika [BMKG] 2019). Another series of significant earthquakes also occurred on August 5 (M w 7.0), August 9 (M w 5.9) and August 19 (M w 6.3 and M w 6.9). Indonesian Agency of Meteorology, Climatology, and Geophysics or BMKG (2019) reported that the hypocenters of these significant earthquakes were mostly concentrated in the north of Lombok Island. The significant earthquakes that occurred several times were indicated as the result of segmentation ruptures formed before the occurrence of a strong earthquake (Wang et al. 2020;Sasmi et al. 2020). Salman et al. (2020) using InSAR data also explained that the significant earthquake sequences came from different hypocenters, and formed a rupture of a thrust fault that tilted to the south. This earthquake sequence might trigger secondary disasters such as landslides (Ferrario 2019). In addition, it is suggested that the 2018 Lombok earthquake might cause a volumetric expansion in the shallow magma plumbing system at Rinjani volcano, which might affect the volcanic activity in the future (Salman et al. 2020).
Information related to the regional stress field around the Lombok Island is less studied, but the strain rate analysis from Bock et al. (2003) shows that the regional strain rate in the Sunda-Banda Arc transition zone is mostly oriented around NE-SW, roughly 45°. On the other hand, the 2018 Lombok earthquake caused fault segmentation with a different strike orientation. The focal mechanism analysis of the regional seismic network (Salman et al. 2020) indicates an east-west relative strike orientation between 80° and 114° of the Lombok earthquake fault system. For this reason, we implemented a crustal anisotropy study to prove whether crustal seismic anisotropy in the Lombok region is more influenced by the regional strain rate or the local fault system associated with the 2018 Lombok earthquake rupture zone. To examine this issue in a greater detail, we utilized 2018 Lombok aftershock data recorded by temporary local seismic stations around this region. The aim of the study is to distinguish between stress and structure-controlled anisotropy at distinct locations located within the 2018 Lombok earthquake area. Thus, this analysis may help to provide better knowledge about crustal anisotropy mechanisms around the region and to substantiate the deformation process related to the 2018 Lombok earthquake.

Seismic network and dataset
The SWS analysis was carried out using waveforms from local events recorded by 16 temporary seismic stations ( Fig. 1) around Lombok Island which were installed during 36 days from 4 August to 9 September 2018 (Sasmi et al. 2020). Most of the events were distributed around LM-09 station which is located near the source of the main earthquake on 19 August. The hypocenter depth of the recorded events varied from 0 to 35 km, and most of them were concentrated at the depth around 25 km (Fig. 1b). The hypocenter locations were initially determined using non-linear method, and then relocated using double-difference method (Sasmi et al. 2020).  Sasmi et al. (2020). The red arrows show compressional strain rate direction from Bock et al. (2003). The colored star symbol shows the hypocentres of significant earthquake. The dashed black line in Fig. 1a represents A-A' hypocentre cross section illustrated in 1b) (Sasmi et al. 2020). The black triangle shows the location of volcanoes. The red triangle depicts Mt. Rinjani. The focal mechanism of the significant earthquake was taken from Global Centroid Moment Tensor (Ekström et al. 2012); b aftershock cross section of A-A' line in Fig. 1a, the width of the window used in cross Sect. 1b is 10 km. The curved dashed line represents the interpretation of the 2018 Lombok earthquake fault from Sasmi et al. (2020). The star symbol depicts the hypocenter of the August 7 (M w 6.9) earthquake; c histogram of the number of earthquake per depth There were several stages of routine data processing that must be prepared before conducting SWS analysis, including selecting waveform data, equalizing sampling frequency, and instrument correction. In the first stage, we selected the data based on these criteria: 1) the data must be recorded by at least three seismic stations; 2) the data have three component seismograms (vertical, E-W, N-S); 3) the waveform data must have clear S wave arrival. This selection process yielded a total of 3112 aftershocks that met the selection criteria for the SWS analysis. The number of S phases from these events was 17,191.
We determined the ray parameters and incidence angles using local 1-D velocity model that had been determined by Sasmi et al. (2020). To avoid the distortion effect of S wave phase at the free surface (Liu et al. 1997), we eliminated the shear waves with incidence angles greater than 40°, a typical cut-off value in SWS analysis (e.g., Araragi et al. 2015;Kanaujia et al. 2019). Nuttli (1961) explained that for a homogeneous half-space with 0.25 Poisson's ratio value, the critical angle is 35° from vertical. In this research, we tolerated this angle to 40° to prevent the rejection of too much useful data. We then applied a 0.2-4 Hz bandpass filter to remove the effect of seismic noise.

Method
To determine the SWS parameters, we initially chose the appropriate waveform window length. The windowing process in the SWS analysis (see Additional file 1: Figure S1) is applied using the short time Fourier transform (STFT) analysis (Nawab and Quaiteri 1988), for the sample area of 0.4 s around the arrival time of the wave phase. The aim of this process is to determine the dominant frequency from each waveform sampling zone (see Additional file 1: Figure S2). From the dominant frequency, the signal period will be obtained which then becomes a reference in the cross-correlation process. The windowing waveform process is one of the important steps in SWS analysis, because the duration length of the waveform affects the cross-correlation results between S fast and S slow wave phases.
An example of a good splitting measurement is depicted in Fig. 2. As shown in Fig. 2, parameters δt and ф were then obtained from the results of cross-correlation or rotation-correlation (RC), to the S fast and S slow phases which have been rotated into the incoming polarization direction and its perpendicular value, for the original filtered waveform. (Bowman and Ando 1987;Vecsey et al. 2008). The cross-correlation process is carried out by varying the values of δt and ф in each of the S wave phase pairs. The increment value of δt used on the crosscorrelation was adjusted to the data sampling interval, which is 100 ms. This process also searched all possible values of ф parameter from 0° to 180° with an interval of 1°. From the cross-correlation results, the correlation coefficient value will be obtained and form a cross-correlation coefficient (CCC) contour pattern. The selected splitting parameters for each event were obtained from the contour coordinates with the highest CCC value, with δt on the y-axis and ф on the x-axis. To maximize the results of SWS parameter determination, there were several criteria implemented: (1) we used a minimum CCC value limit of 0.55. We take this value to avoid too much data being wasted due to the tight selection of CCC minimum threshold. Only values above a certain threshold are used in further processing; (2) the particle motion of the corrected waveforms has a linear pattern after correction; (3) a good standard deviation for the ф parameter is usually less than 30°, and for the δt is 0.02 s (e.g., Hurd & Bohnhoff 2012).

Delay time and fast polarization analysis
There are several factors that affect the variation of the SWS polarization direction value, such as the dimensions and orientation of the subsurface anisotropic material, the presence of anisotropic mineral layers, the intensity of fractures in the medium, the presence of rock layers, etc. (e.g., Crampin et al. 1986;Shelley et al. 2014). In the earth's crust, the direction of ф and δt is generally very dependent on the presence of fractures, sedimentary layers, and the direction of the maximum horizontal stress. In tectonic active areas, the direction of S fast tends to be parallel to the horizontal direction of the maximum stress (Crampin and Chastin 2003). In addition, paleostress conditions can also contribute to the direction of ф (Blenkinsop 1990). Kleinrock and Humphris (1996) stated that this condition occurs because the stress field tends to cause joint alignments, rock pores, and microcracks which can increase the degree of seismic anisotropy. The value of the δt is strongly controlled by the degree of anisotropy of the medium and also the intensity of the fracture. The δt will be smaller if the fracture density in the medium is low.
The average ф and δt values obtained from this study are summarized in Table 1. The SWS circular statistical analysis to determine the average and standard deviation of ф was performed following Mardia (1972) and Davis (1986). Based on the rose diagram ( Fig. 3) and Table 1, we observe that the fast polarizations at most observation stations tend to vary in some different areas. Commonly, in the earth's crust, the ф at an observation station is relatively parallel to the strike of the local fault (e.g., Eken et al. 2011;Zhang et al. 2018;Audet 2014;Boness and Zoback 2006). This condition can be seen at several stations adjacent to the fault in this study, such as     (2020) estimated the fault orientation of the Lombok earthquake rupture area using the waveform inversion of the regional seismic network, and obtained the strike direction range of ф similar to the ф of our result. In that study it was explained that the July 28 earthquake had a strike plane of the N114°E fault; the strike of the August 5 earthquake fault was N80°E; August 19 earthquake Mw 6.3 and 6.9 had strikes of N100°E and N95°E, respectively. Thus, it can be concluded that the SWS parameters at several stations distributed in North Lombok are related to the presence of the 2018 Lombok earthquake faults system. Analyzing the character of the delay time with respect to depth can provide general information about the position of the anisotropic medium (e.g., Hurd & Bohnhoff 2012;Syuhada et al. 2017;Hiramatsu and Idaka, 2015). In Fig. 1c, we can see that the earthquake hypocenters are mainly distributed to the depth of 15-20 km. Figure 4a shows δt variation for every depth increment. This figure depicts that the higher δt is mainly distributed at the first 15 km depth, with the highest value of δt up to 0.38 s. This pattern of the δt also shows that δt trend increases at depths 15-20 km, then apparently decreases after that depth as the data amount reduces. The condition suggests that the possible dominant source of anisotropy is located in the upper crust zone, approximately 15-20 km. Syuhada et al. (2017) also found through SWS analysis using the regional seismic station in the Sunda-Banda arc transition zone that the anisotropy in that region is confined down to 15-25 km depth. Meanwhile, Kaneshima (1990) in his study about crustal anisotropy in Japan also reported higher δt observed for events with focal depth up to 15 km, suggesting that crustal anisotropy existed in the upper 15 to 25 km depth of the crust. Meanwhile, the decreasing trend of δt in the depth of > 20 km may indicate a decrease in the anisotropy structure, for example due to the closing of the fault zone or the thinning of the mineral layer zone (Kaneshima 1990). On the other hand, variation of ф with respect to depth (Fig. 4b) does not show any particular observable trend. An explanation of this condition in the next chapter will be described in Discussion chapter.

Spatial averaging of SWS fast polarization direction
To substantiate and analyze the spatial variations of fast polarization direction, we adapted the shear wave splitting spatial averages (TESSA) program algorithm developed by Johnson et al. (2011);Johnson and Savage (2012). In this case, the ф values were determined to every spatial grid block passed by a certain number of rays, and weighted inversely to the square of the distance from the station. The minimum block size assigned was 10 km containing 10-80 rays in each grid block. If the standard deviation was greater than 30° or the raypath density was too low (< 10), the result of ф spatial average would not be plotted in the block. The mean fast direction ф was then calculated using circular statistics (Mardia 1972;Davis 1986) (see Additional file 1: Figure S2).
We analyzed the complexity of the spatially averaged ф, and the result is shown by Fig. 5. Through the analysis of the spatial average of ф, we will further analyze the causes of variations in the direction of fast polarization in this study, whether it is controlled by regional stress conditions or local geological structures (e.g., Hiramatsu and Iidaka 2015).
SWS spatial analysis needs to be done by separating the data per depth based on the area of the crust or under the crust to make it easier to analyze the anisotropic structure causing the splitting (e.g., Syuhada et al. 2017;Hiramatsu et al. 2010). Based on information from Curray et al. (1977) and Yang et al. (2020), the depth of the crust in the Lombok area ranges from 20-30 km. Considering the reference, we constrained the depth up to 35 km to perform the spatial average of SWS.
In Lombok strike slip fault and Sumbawa strike slip fault (Fig. 5), the ф values are generally parallel to the strike of the two faults, which is around N 30°E (Irsyam et al. 2017). Around Sumbawa Island particularly in the oceanic region, the spatial average of fast directions is also in the NE-SW direction. This condition suggests that crustal anisotropy around this area may still be controlled by the opening of microcracks related to the local tectonic stress (Syuhada et al. 2017;Bock et al. 2003), but it does not rule out the possibility that local structures such as microcracks can still contribute to the distribution of the ф value. Meanwhile towards the northern region of Lombok Island, the spatially averaged fast polarizations are predominantly around W-E. This condition most likely occurs due to the lineament of the Flores Back Arc Thrust, or the fault segmentation system caused by the 2018 Lombok earthquake, which shows a similar lineament (Salman et al. 2020). This resemblant pattern can be observed in the seismograph stations located to the north of Lombok Island (LM 03, LM 04, LM 05, LM 06). The further explanation of our findings is described in Discussion chapter.

Temporal variation analysis
Analysis of temporal variation in splitting parameters is a recommended method for predicting the occurrence of earthquakes (e.g., Aster and Shearer 1991;Syuhada et al. 2017) and volcanic eruptions (e.g., Miller and Savage 2001;Gerst and Savage 2004;Gerst 2003). The microcrack geometry changing prior to an earthquake or volcanic eruption can affect the change in splitting parameters.  reported that the temporal variations in δt have been observed after the occurrence of large earthquakes showing the systematic pattern between the event magnitudes and delay time. Saiga et al. (2003) also reported that there is a relation between the changes in SWS parameters and the expanding crack density due to the coseismic stress changes after the Aichi-Ken Tobu earthquake, Japan. Furthermore, temporal variations of ф have been utilized as a possible tool to analyze the stress change beneath Mt. Ruapehu, New Zealand, prior to the volcanic eruption (Gerst and Savage 2004). This previous study suggested that the changes in splitting parameters around this region are associated with the pressure changes originating from the magma chamber. Temporal variations in SWS parameters can also be employed to detect the fluid-saturated microcracks zone , and the alteration of microcrack geometry related to the intermediate earthquakes .
We utilized data from 36 days of earthquake recording time from August 4-September 9, 2018, to perform the temporal analysis of SWS parameter. We display two examples of stations that probably exhibit temporal variations in SWS parameters at LM05 and LM06 stations. Figure 6 shows the temporal variation in time delay at station LM05 and LM06. We applied 10-window-moving averages of splitting delay time. We observe interesting patterns of the delay-time moving-average that are probably correlated to the magnitude variation of the events. The two stations indicate correlation between temporal variations in δt and the earthquake magnitudes. In general, the pattern of δt seems to increase slightly following the significant earthquake on August 9th. The increasing trend is also clearly seen after the August 19 earthquake, given that on August 19, 2018 there were two fairly large earthquakes (M w 6.3 and M w 6.9).

Discussion
Tectonically, the Lombok region is located in the transition zone of the Sunda-Banda arc, so that there are many fault systems in the area. It is possible for microcracks to be distributed from various fault sources (Wang et al. 2020). However, in this study, the earthquake sources used in the SWS data processing were mostly from the aftershock zone of the Lombok earthquake. This indicates that the SWS parameters observed in this study most likely reflect the anisotropic character due to the Lombok Earthquake fault system. In addition, one of the assumptions of the SWS method is that the ф reflects the last interval of anisotropy encountered along the raypath, and the δt parameter represents the amount of anisotropy accumulated along the entire ray path (Yardley and Crampin 1991). Therefore, it can be assumed that in this study, the local structure around Lombok Island may contribute significantly to the SWS parameters distribution in the area's local seismic network.
Previous research in the Lombok Island area and the Sunda-Banda arc transition zone showed that the SWS parameters in the region were still influenced by the stress-induced anisotropy factor due to plate movement (Syuhada et al. 2017), where the direction of the regional strain rate axis in the Lombok area is generally oriented in NE-SW direction (Bock et al. 2003), and the previous station which is located in the north of Lombok showed fast polarization average of NNE-SSW (Syuhada et al. 2017). Using local station network data of the 2018 earthquake aftershock, we can understand that in north Lombok, the seismic stations tend to head in the E-W direction with an average fast polarization value of 85.48-93.03°. This may occur due to the development of the fault system due to the 2018 Lombok earthquake, which tends to have an E-W direction, so it affects the change in the direction of S fast polarization. Thus, it is possible that the mechanism of SWS parameters in the Lombok area underwent a significant change after the 2018 Lombok earthquake, as inferred from our result of spatial averaging of ф.
Through the SWS analysis, we become more confident that there is a large surface deformation after the 2018 Lombok earthquake in the northern area of Lombok Island, which affects the anisotropic character of the SWS parameters. The information from our study is expected to be useful for future studies, such as focal mechanism and stress analyses during the 2018 Lombok earthquake.
The fast polarization spatial averaging analysis (Fig. 5) can help us understand the relationship between fast polarization parameters and local structures. We can observe the similarity of the pattern between the direction of the local structure's lineaments and the direction of the ф parameter at adjacent stations. In the central area of Lombok Island adjacent to Mount Rinjani, the pattern of ф direction is slightly various from NW-SE to NE-SW. Johnson et al. (2011) and Shelley et al. (2014) explained that in a volcanic area, the stress field tends to be perturbed by magmatic intrusion and migration in the lithosphere, so that it can affect the distribution and the direction of crack within the crust.
Thus, we estimate that the ф variation in this area probably reflects the anisotropic character due to the presence of volcanic activity in Mount Rinjani. In addition, the complexity of the anisotropy sources in this region may cause the Rosette Diagram to spread in various directions (Johnson 2013). This condition is also illustrated in seismic station LM07 where the rosette diagram of ф shows a distribution over a wide range of values.
The time delays observed between the slow and fast shear waves give information about the density of cracks in the medium. The relation between δt and depth shows interesting pattern of possible anisotropy structure Fig. 6 a Temporal variation in delay time at station LM05; b temporal variation in delay time at station LM06; c variation of earthquake magnitude with respect to recording time. The grey dots in Fig. 6a and 6b show the delay time variation before moving-average is applied. The green dots in Fig. 6a and 6b show the delay time variation after 10-window moving-average. Blue dots in Fig. 6c show the magnitude variation before the moving-average process is applied. Orange dots in Fig. 6c show the magnitude variation after generate the 10-window moving-average. Blue oval symbols depict the trend correlation with respect to August 9 (M w 5.9) earthquake, and yellow oval symbols show the correlation towards August 19 earthquakes magnitude (Fig. 4a). There is an increase trend of δt in the depth of < 15-20 km that probably related to the influence of the crack zone formed during the 2018 Lombok significant earthquake (Sasmi et al. 2020;Afif et al. 2021;Priyono et al. 2021;Wang et al. 2020;Salman et al. 2020). Sasmi et al. (2020) explained that the 2018 Lombok earthquake aftershocks were mostly distributed to the depth of < 20 km with an average hypocenter slope of ~ 30° and thought to be related to the rupture zone of 2018 Lombok earthquake. These significant events of earthquake were suspected as a result of the oceanic crust subducted to the north of Lombok Islands and caused fault segmentation (Sasmi et al. 2020;Afif et al. 2021;Priyono et al. 2021;Wang et al. 2020;Salman et al. 2020;Yang et al. 2020;Lythgoe et al. 2021). The results from these studies confirm the shallow anisotropic fabrics may control the main possible sources of anisotropy in this area (e.g., Syuhada et al. 2017).
Other than that, the delay time trend also rises specifically at the depth of < 15 km (Fig. 4). We interpret that this condition indicates the continuity of the 2018 Lombok earthquake fault from the hypocenter depth to a shallower surface. Crampin et al. (1984); Crampin and Booth (1985) explained that the presence of pore-fluid can cause the fracture critically of rocks and generate anisotropy of 0.5-10% in the upper crust. We suggest that a fluid inclusion, either from fluid-filled cracks from Mt. Rinjani volcanic system (Afif et al. 2021;Sundhoro et al. 2000), or fluid related to the subduction system of oceanic crust in the north of Lombok (Afif et al. 2021;Priyono et al. 2021;Yang et al. 2020) contributes to the crustal anisotropy beneath the Lombok Island. The presence of fluid in the microcrack systems can preclude the microcracks to close producing higher degree of anisotropy. This interpretation can be supported by the high ratio of Vp/Vs (Afif et al. 2021) and the value of the seismic attenuation parameter (Qp) which is quite low beneath Mount Rinjani (Priyono et al. 2021). The high temperature of the area is also shown by the attenuation tomography result (Priyono et al. 2021). This condition of high ratio of Vp/Vs and low Qp parameter beneath a volcanic area indicates that this area may contain fluids, either high-pressure hydrothermal water or melt (e.g., Ramdhan et al. 2019;Nakajima and Hasegawa 2006). The fracture system of depth < 15 km is interpreted as a pathway of fluid to rise to the surface. The remaining fluids in the study area are suspected to be related to the subsurface geothermal system, considering that in the Lombok area there is Mount Rinjani as one of the controllers of the local geological setting. Sundhoro et al. (2000) described the geothermal manifestations and surface alteration in the north of Mount Rinjani, in the form of hot springs at Lake Segara Anakan, Kukusan, Kalak, and Sebau. The alleged presence of fluid in the presumably permeable zone at a depth of < 15 km may be related to the flow of hydrothermal fluid supply from the subsurface to the surface so as to form geothermal manifestations around Mount Rinjani.
Delay time variation generally decreases with increasing depth. In our case, the variation of the delay time value also begins to decrease after a depth of > 20 km. The decrease in delay time with depth occurs because as the depth increases, the fracture will close due to an increase in overburden pressure. In certain cases, an increase in delay time at large depths can occur due to the presence of fracture filling fluid at that depth. The fluid can increase the pore pressure and prevent the fracture from closing again due to high overburden pressure (e.g., Syuhada et al. 2017;Hiramatsu and Iidaka 2015). However, we still need more findings to prove our interpretation. This provisional result can be a reference for performing SWS delay time tomography to prove this hypothesis.
Meanwhile, we do not discuss in more detail about the pattern of the ф parameter with respect to depth because we did not find a specific trend in the variation of with depth in this study (Fig. 4b). This condition of scattered ф with respect to depth is common in the earth's crust area (e.g., Syuhada et al. 2017;Martin et al. 2014), because the variation in the ф value is still very much controlled by the relatively complex anisotropy of the earth's crust. For example, this is due to cracks that have various directions, anisotropic sediment and mineral layer structures, microcracks due to stress-induced anisotropy, or a combination of these causes.
As mentioned before, the magnitude of the ф parameter is related to the direction of the Sfast wave polarization, which also reflects the condition of the structure of the seismic anisotropy medium, including microcracks, fault structures, and anisotropic minerals (Savage et al. 2016). In general, the value of the ф parameter in an area is most likely sub-parallel with the strike direction of the dominating anisotropic structure in that area.
In the Lombok region, the direction of the strain rate tends to be NNE-SSW (Syuhada et al. 2017;Bock et al. 2003). Meanwhile, Irsyam et al. (2017) report that the Lombok and Sumbawa strike-slip fault structures tend to be oriented NE to SW. Salman et al. (2020) explained that there is a fault structure in the north of Lombok Island with a strike around N 88-93° E, which might be related to the Lombok earthquake fault system. The similar results are also observed in our splitting results through the mean ф orientation plots (Fig. 3) and also the ф spatial distribution map (Fig. 5). However, the plot of ф versus delay time in Fig. 4b does not reflect the dominance of a particular ф value trend. We suspect the absence of the trend is because the ф versus delay time graphic reflects the accumulation of all ф values in the entire area and throughout the depth, so the range of values is wide. Thus, the plot of the trend of delay time over time tends not to be dominated by certain structural directions down to the depth of 35 km.
Scatter conditions in the delay time versus depth plot have also been found in several studies in the crustal region, such as Syuhada et al. (2017) and Hiramatsu et al. (2010). Suppose the range of ф values recorded in that area is wide and spread out enough. In that case, it can be assumed that the area has a complex subsurface anisotropic structure, which also can be observed in our study, as shown in Fig. 4b. Therefore, we suggest that the absence of the trend of the ф parameter with respect to depth is due to the complexity of the anisotropic structure of the study area due to the Lombok earthquake fault.
As we have implemented and explained before, spatial analysis of the SWS parameter can help to provide general information about the location of the anisotropic zone. Meanwhile, we can also understand when the expansion of the anisotropic zone occurs that might be related to the increase of fracture density or crack aspect ratio through a temporal analysis of the SWS parameter. Liu et al. (1997) used temporal analysis of SWS to determine the relationship between changes in delay time trend with increased seismic activity around Parkfield, California. Miller and Savage (2001) suggested SWS temporal analysis as a recommended method to analyze the period of microcrack formation that can be utilized to predict an upcoming volcanic eruption.
Temporal analysis of the delay time parameter shows that the delay time trend varies throughout the time for both of the two stations (Fig. 6). It can be observed that there is a potential correlation between the movingaverage of delay time and the earthquake magnitude. The general trend of delay time is seen to increase following the rising trend of earthquake magnitude over time. The enhancement of the trend appears to be quite significant around the occurrence of the Mw 6.3 and Mw 6.9 earthquakes on August 19, 2018.  suggested that the increase in delay time may indicate a crack aspect ratio change related to the stress change that occurs after a strong earthquake event.
To further analyze this condition, we mapped the spatial variation of the hypocenter with respect to time, as shown in Fig. 7. The temporal variation of the hypocenter initially showed a more focused distribution of the earthquake hypocenter in the north and west of Lombok Island at the period 4-17 August 2018 (Sasmi et al. 2020). At around August 18 and 19, a new aftershock cluster was formed around the east to northeast of Lombok Island, indicating the presence of the large deformation formed after the August 19 earthquake around the east to the northeast part of Lombok Island (Sasmi et al. 2020;Supendi et al. 2020;Wang et al. 2020). Generally, the increasing pattern in delay time variation with respect to the time (Fig. 6) may be related to an increase in the crack density pertained to the stress change subsequent to the significant earthquake  of the 2018 Lombok earthquake. Wang et al. (2020) explained that there is a stress increase surrounding the east of the August 19 seismogenic fault, as shown by the result of the Coulomb stress analysis. The stress was thought to be accumulated before the August 19 earthquake and then released, resulting in slip movement generating M w 6.3 and M w 6.9 earthquakes (Wang et al. 2020;Salman et al. 2020).
On the other hand, the condition of the data is not yet sufficient to properly describe the temporal variation around August 7th. This was due to the incomplete seismographic network installed at the location on that date. Although we are able to examine the variation of δt correlated with the event magnitudes, the temporal trend of ф is not obviously seen in any station. Similar findings are also reported by Syuhada et al. (2017), which analyzed the temporal variation of SWS parameters in adjacent areas using regional seismographic network data. Yang et al. (2011) also did not find any clear trend of ф in the SWS temporal variation of southern California. As explained by Nur (1971), delay time generally depends on crack density, crack aspect ratio, and anisotropic path length, while the fast polarization direction depends on the direction of the horizontal compressional strain rate and fault strike. Several studies (Crampin and Chastin 2003;Zatsepin and Crampin 1997) show that before a significant earthquake occurs, the stress is accumulated until it reaches a critical level for the rock to fracture and lose its shear strength. This is the reason why earthquakes can be monitored through the changes in local stress concentrations. In SWS analysis, this change causes an increase in crack density and crack aspect ratio and can be identified by an increase in the average delay time over time. In this case, the character of the temporal magnitude change shows a stronger relationship with the delay time rather than with the fast polarization direction. However, the provisional assumption in our result requires more data and further research to prove the actual occurrence of the fault continuity to the surface.
One of the difficulties in measuring the SWS phenomenon on the small-magnitude aftershocks is the complicated nature of the waveform, which causes the spatial and temporal δt variations to scatter, which generally ranges from ± 80% of the mean (Crampin and Gao 2006).
We can observe this condition in Fig. 6a and Fig. 6b, where the scatter on the parameter δt is very visible at stations LM05 and LM06 before moving averages are carried out on the data to reduce the bias effect due to the spreading data. Volti and Crampin (2003) explained that the cause of this scatter may be due to several things: anisotropic variations in direction, errors in determining the location of the earthquake hypocenter, errors in determining the arrival time of the S fast and S slow waves, complex geological conditions, and complicated crack geometry. We have carried out the process of determining the arrival time of waves and relocating the hypocenter carefully (Sasmi et al. 2020). Then we suspect that the complex anisotropic conditions may contribute to the delay time variation in our results.
However, distinguishing between the influence of stress-and structure-induced anisotropy is quite tricky. Structural-induced anisotropy is mostly caused by the presence of local fault zones or foliation of anisotropic medium, so that the fast polarization will be parallel to the strike of the anisotropic structure. It is said to be tricky if there is an anisotropic structure that has a strike that is close to the compression direction of the strain rate area, so it will be difficult to distinguish the cause of the anisotropy (e.g., Johnson et al. 2011). Since the local structure may contribute to anisotropy significantly, the splitting parameters obtained at stations located to the structure indicates the fast direction to be parallel to the plane of the structure or mineral alignment and can be different to the tectonic stress orientation. As shown in Fig. 5, most likely the SWS parameters are strongly controlled by the local structure rather than the regional compressional strain.
Otherwise, to observe the stress-induced anisotropy, we may need more stations that are located further from any active structure (e.g., Li and Peng 2017). However, in this study, we do not have data from seismic stations outside the Lombok fault zone, so our  Sasmi et al. (2020) in each time sequence shown in the white title. The inverted triangles show the location of seismometers. Red triangle shows Mt. Rinjani. The yellow diamond depicts the hypocenter of July 29 (M w 6.4) earthquake (BMKG, 2019). The yellow stars represent the significant earthquake of August 5 (M w 7.0), August 9 (M w 5.9), and 19 (M w 6.3 and 6.9) determined by Sasmi et al. (2020) interpretation is limited to seismic stations located within the interest area of 2018 Lombok earthquake.

Conclusion
SWS parameter analysis and spatial averaging of fast azimuths have been carried out around the hypocenter area of the 2018 Lombok earthquake aftershock using 16 local seismographic network data. The results of the SWS measurements show that the ф distribution implies consistency with the presence of local fault structures, such as the Lombok Shear Fault, Sumbawa Shear Fault, and the fault area due to the Lombok earthquake 2018. We also observed a relatively E-W trend of ф direction in the north of Lombok Island, which is consistent with the strike direction of fault formed after the 2018 Lombok Earthquake, and Flores Back Arc Thrust. The NNE fast azimuths observed in the west and east part of study area may also be associated with the presence of Lombok Strait Strike Slip Fault and Sumbawa Strait Strike Slip Fault. We suggest that the anisotropy sources beneath the 2018 Lombok Earthquake rupture zone are mainly controlled by the presence of local structure rather than the compressional stress setting of the area, so that the SWS parameter in the 2018 Lombok Earthquake aftershock zone can be categorized as structure-induced anisotropy. The anisotropy condition around the upper crustal area with a depth of < 15 km is allegedly controlled by fluid inclusion, which is suspected either from Mt. Rinjani magmatic activity or crack-filled fluid of Flores Oceanic Crust subduction system near the main source of the 2018 Lombok Earthquake. The δt versus depth graphic shows that the anisotropy structure is likely found in 15-20 km depth and consistent to Lombok earthquake source depth. The temporal variation in δt observed in this study may be correlated with the stress changes following the increasing trend of event magnitudes. Therefore, the anisotropy mechanism in this area is strongly presumed to be related to the structure-controlled anisotropy.
Additional file 1: Figure S1. Example of window selection using STFT for waveform cross-correlation process. a). An example of S fast waveform which is sampled with 0.4 s window around the arrival time area. b). The spectrum of S2.a). waveform with a dominant frequency of 4.7 Hz; c). An example of S slow waveform which is sampled with 0.4 s window around the arrival time area; d). The spectrum of S2.b). waveform with a dominant frequency of 4.7 Hz. Figure S2. The map depicts the grid blocks (grey lines) and rays (blue lines) used in the spatial averaging of ф from the temporary deployment data. The grid increments used in the inversion is 10 x 10 km.