Spatio-Temporal Distribution of Ground Deformation Due to 2018 Lombok Earthquake Series

Lombok Island in Indonesia was hit by four major earthquakes (6.4 Mw to 7 Mw) and by at least 818 earthquakes between 29 July and 31 August 2018. The aims of this study are to measure ground deformation due to the 2018 Lombok earthquake series and to map its spatio-temporal distribution. The application of DinSAR was performed to produce an interferogram and deformation map. Time series Sentinel-1 satellite imageries were used as master and slave for each of these four major earthquakes. The spatio-temporal distribution of the ground deformation was analyzed using a zonal statistics algorithm in GIS. It focused on the overlapping area between the raster layer of the deformation map and the polygon layer of six observation sites (Mataram City, Pamenang, Tampes, Sukadana, Sembalun, and Belanting). The results showed that the deformation includes uplift and subsidence. The first 6.4 Mw foreshock hitting on 29 July 2018 produces a minimum uplift effect on the island. The 7.0 Mw mainshock on 5 August 2018 causes extreme uplift at the northern shore. The 6.2 Mw Aftershock on 9 August 2018 generates subsidence throughout the study area. The final earthquake of 6.9 Mw on 19 August 2018 initiates massive uplift in the study area and extreme uplift at the northeastern shore. The highest uplift reaches 0.713 m at the northern shore, while the deepest subsidence is measured −0.338 m at the northwestern shore. Dominant deformation on the northern area of Lombok Island indicates movement of Back Arc Trust in the north of the island. The output of this study would be valuable to local authorities to evaluate existing earthquake’s impacts and to design mitigation strategies to face earthquake-induced ground displacement.


Introduction
DInSAR (Differential Interferometric Synthetic Aperture Radar) is a satellite-based approach in Geographic Information Science that allows measuring of deformation on earthquake-impacted areas [1][2][3][4] in the last decades [5,6]. Integrating DinSAR and GIS would improve the studies related to topographic changes using conventional technique [7][8][9][10] and direct measurement [11][12][13][14][15]. DInSAR can be created from active remote sensing that transmits radio waves with vertical and horizontal polarization to the earth's surface and then receives it back to form radar images. The main principle of DInSAR is to calculate the interference between two radar images and measure the displacement value between them. Most of the studies extensively used the DinSAR technique to analyze a single earthquake event, such as in Chile [16,17], New Zealand [18], Iran-Iraq [19][20][21], US [22], Cyprus [23], Azerbaijan [24], Italy [25], and Mexico [26]. Although, the occurrence of earthquake series should be studied as well.
Sentinel 1 allows to measure the displacement due to earthquake [2,3,5]. It was launched in 2014 as a continuation of the European Space Agencies' previous SAR satellite channels ERS and ENVISAT. Its mission is dedicated to providing high-resolution remote sensing image data to support observation of environmental changes on earth's surface to ensure sustainable economic growth. Sentinel 1 operates C-band sensor (5.405 GHz) installed on two satellite constellations (Sentinel-1A and Sentinel-1B) with repeat cycles every six days. Sentinel 1 is working on both single polarization (HH and VV) and double polarization (VV + VH and HH + HV). This outstanding specification is then fit into the data requirement in order to monitor ground displacement due to earthquake series.
Lombok Island in Indonesia experienced a massive earthquake series in 2018 ( Figure 1). This island is between Back Arc Thrust ( Figure 2) in the north and active subduction zone of the Java trench in the south due to the convergence movement of the Eurasian and Indo-Australian plates. Four tertiary formations from early Miocene to late Miocene is located in the southern part of the island where most of the lineaments can be found. These geological formations include Kawangan formation (alternating quartz sandstone, claystone, and breccia), Pengulung formation (breccia, lava, tuff with lenses of limestone), Ekas Formation (Limestone/calcarenite), as well as Dacite and Basalt intrusive rocks [27]. Four younger formations from Pleistocene are located in the middle of the island. They consist of Kalipalung formation (alternating calcareous breccia and lava), Selayar Formation (tuffaceous sandstone, tuffaceous claystone with thin carbon internalations), Kalibabak formation (breccia and lava) and Lekopiko formation (pumiceous tuff, lahar breccia and lava). The youngest formations in the northern part of the island are from Holocene. At low altitude, there is Alluvium (pebble, granule, sand, clay, and coral fragments), while at high there is undifferentiated volcanic rocks originated from Mount Pusuk-Nangi and from Mount Rinjani. The latter is the present name of Samalas volcano which experienced caldera-forming eruption in 1257 AD, producing 40 km 3 volcanic materials including deposits in both the Arctic and Antarctic ice cores [28]. This volcano is located next to Back Arc Thrust in the north of the island. An interesting finding related to the location of the Back Arc Thrust is that its location is actually much closer to the shore than was expected from the previous geological map [1] (Figure 2). between them. Most of the studies extensively used the DinSAR technique to analyze a single earthquake event, such as in Chile [16,17], New Zealand [18], Iran-Iraq [19][20][21], US [22], Cyprus [23], Azerbaijan [24], Italy [25], and Mexico [26]. Although, the occurrence of earthquake series should be studied as well. Sentinel 1 allows to measure the displacement due to earthquake [2,3,5]. It was launched in 2014 as a continuation of the European Space Agencies' previous SAR satellite channels ERS and ENVISAT. Its mission is dedicated to providing high-resolution remote sensing image data to support observation of environmental changes on earth's surface to ensure sustainable economic growth. Sentinel 1 operates C-band sensor (5.405 GHz) installed on two satellite constellations (Sentinel-1A and Sentinel-1B) with repeat cycles every six days. Sentinel 1 is working on both single polarization (HH and VV) and double polarization (VV + VH and HH + HV). This outstanding specification is then fit into the data requirement in order to monitor ground displacement due to earthquake series.
Lombok Island in Indonesia experienced a massive earthquake series in 2018 ( Figure  1). This island is between Back Arc Thrust ( Figure 2) in the north and active subduction zone of the Java trench in the south due to the convergence movement of the Eurasian and Indo-Australian plates. Four tertiary formations from early Miocene to late Miocene is located in the southern part of the island where most of the lineaments can be found. These geological formations include Kawangan formation (alternating quartz sandstone, claystone, and breccia), Pengulung formation (breccia, lava, tuff with lenses of limestone), Ekas Formation (Limestone/calcarenite), as well as Dacite and Basalt intrusive rocks [27]. Four younger formations from Pleistocene are located in the middle of the island. They consist of Kalipalung formation (alternating calcareous breccia and lava), Selayar Formation (tuffaceous sandstone, tuffaceous claystone with thin carbon internalations), Kalibabak formation (breccia and lava) and Lekopiko formation (pumiceous tuff, lahar breccia and lava). The youngest formations in the northern part of the island are from Holocene. At low altitude, there is Alluvium (pebble, granule, sand, clay, and coral fragments), while at high there is undifferentiated volcanic rocks originated from Mount Pusuk-Nangi and from Mount Rinjani. The latter is the present name of Samalas volcano which experienced caldera-forming eruption in 1257 AD, producing 40 km 3 volcanic materials including deposits in both the Arctic and Antarctic ice cores [28]. This volcano is located next to Back Arc Thrust in the north of the island. An interesting finding related to the location of the Back Arc Thrust is that its location is actually much closer to the shore than was expected from the previous geological map [1] (Figure 2).    [27]). Back Arc Thrust is located next to the northern shore [1]. Photos of building structural damage were taken at alluvial area of Mataram City.
This geologic setting leads to continuous exposure of this island to earthquakes. Recent shallow and great magnitude of earthquake series occurred in 2018 next to the Back Arc Thrust. The first 6.4 Mw foreshock hit on 29 July 2018. It was followed by a 7.0 Mw mainshock on 5 August 2018, then a 6.2 Mw aftershock on 9 August 2018, and finished by an earthquake of 6.9 Mw on 19 August 2018 ( Figure 1). In between and after these four days of strong motion occurrence, there were at least 818 minor earthquakes until 31 August 2018 at 18:00 Western Indonesian Time. This led to 560 deaths, 396,032 evacuated people, and 83,392 destructed houses (Figure 2). Such a great number of destructed houses due to the high magnitude of consecutive earthquakes in a short period of time would be caused by the abrupt change in local topography. This is well-known as ground deformation which may involves uplift and/or subsidence. The aims of this study are to measure ground deformation due to the 2018 Lombok earthquake series and to map its spatiotemporal distribution. This would be important to local authorities to evaluate the earthquake's impacts which still exist now and mitigate future disasters related to earthquakeinduced ground displacement.

Materials and Methods
This study was carried out on Lombok Island, Indonesia using Sentinel 1 satellite imageries. Indeed, the application of Sentinel 1 has been used worldwide for deformation analysis mainly due to subsidence in metropolitan areas [29][30][31][32][33][34], volcanic activity [35][36][37][38][39], landslides [29,[40][41][42][43][44], and single-event earthquakes [17][18][19][20][21][22]. In order to study the earthquake series in Lombok, the methodology of this research focused on a DInSAR creation in SNAP using Two Pass Interferometry method based on 2 radar images with different recording times ( Figure 3). One image acted as the master image and another image as the slave image (Table 1). These two images were paired to form an interferogram to produce a deformation map with uplift or subsidence values. Data selection depended on the mode and intensity of the images which was related to image polarization [2,5]. The first stage of radar processing was TOPSAR Split, selecting bursts from Sentinel 1 radar imagery in the study area. There are three types of modes in each image obtained from the TOPSAR Figure 2. Geology map of Lombok Island (adapted from [27]). Back Arc Thrust is located next to the northern shore [1]. Photos of building structural damage were taken at alluvial area of Mataram City.
This geologic setting leads to continuous exposure of this island to earthquakes. Recent shallow and great magnitude of earthquake series occurred in 2018 next to the Back Arc Thrust. The first 6.4 Mw foreshock hit on 29 July 2018. It was followed by a 7.0 Mw mainshock on 5 August 2018, then a 6.2 Mw aftershock on 9 August 2018, and finished by an earthquake of 6.9 Mw on 19 August 2018 ( Figure 1). In between and after these four days of strong motion occurrence, there were at least 818 minor earthquakes until 31 August 2018 at 18:00 Western Indonesian Time. This led to 560 deaths, 396,032 evacuated people, and 83,392 destructed houses ( Figure 2). Such a great number of destructed houses due to the high magnitude of consecutive earthquakes in a short period of time would be caused by the abrupt change in local topography. This is well-known as ground deformation which may involves uplift and/or subsidence. The aims of this study are to measure ground deformation due to the 2018 Lombok earthquake series and to map its spatio-temporal distribution. This would be important to local authorities to evaluate the earthquake's impacts which still exist now and mitigate future disasters related to earthquake-induced ground displacement.

Materials and Methods
This study was carried out on Lombok Island, Indonesia using Sentinel 1 satellite imageries. Indeed, the application of Sentinel 1 has been used worldwide for deformation analysis mainly due to subsidence in metropolitan areas [29][30][31][32][33][34], volcanic activity [35][36][37][38][39], landslides [29,[40][41][42][43][44], and single-event earthquakes [17][18][19][20][21][22]. In order to study the earthquake series in Lombok, the methodology of this research focused on a DInSAR creation in SNAP using Two Pass Interferometry method based on 2 radar images with different recording times ( Figure 3). One image acted as the master image and another image as the slave image (Table 1). These two images were paired to form an interferogram to produce a deformation map with uplift or subsidence values. Data selection depended on the mode and intensity of the images which was related to image polarization [2,5]. The first stage of radar processing was TOPSAR Split, selecting bursts from Sentinel 1 radar imagery in the study area. There are three types of modes in each image obtained from the TOPSAR technique, i.e., IW1, IW2, and IW3. The choice of mode type depended on the research area. All modes in each image generally had different observation areas. In this study, we chose three bursts in IW3 to observe the research area. Perpendicular baseline affects the signal Remote Sens. 2021, 13, 2222 4 of 14 phase for deformation detection. Too short of a baseline reduced the sensitivity of signal phase, while too long of a baseline produced noise resulting in phase errors [45]. technique, i.e., IW1, IW2, and IW3. The choice of mode type depended on the research area. All modes in each image generally had different observation areas. In this study, we chose three bursts in IW3 to observe the research area. Perpendicular baseline affects the signal phase for deformation detection. Too short of a baseline reduced the sensitivity of signal phase, while too long of a baseline produced noise resulting in phase errors [45]. The intensity of the signal was checked on one type of twin vertical polarization (VV). VV polarization was chosen because it proposed the most ideal observation of the deformation of the earth's surface since the angle of signal transmission and reception focused on vertical observations. Thus, the intensity of VV polarization was greater than that of VH and HV polarization. We also applied Orbit File in order to conduct corrections to each image. Sentinel Precise Orbits contained information about satellite positions during SAR acquisition. Precise Orbit Ephemerides (POE) covered approximately 28 hours and contained orbit state vectors at fixed time steps with 10 second intervals. This orbit file was generated per day and sent within 20 days of data acquisition. It was important to unveil the information on the data type, acquisition date, position, and orbit of the images during data acquisition.   The intensity of the signal was checked on one type of twin vertical polarization (VV). VV polarization was chosen because it proposed the most ideal observation of the deformation of the earth's surface since the angle of signal transmission and reception focused on vertical observations. Thus, the intensity of VV polarization was greater than that of VH and HV polarization. We also applied Orbit File in order to conduct corrections to each image. Sentinel Precise Orbits contained information about satellite positions during SAR acquisition. Precise Orbit Ephemerides (POE) covered approximately 28 h and contained orbit state vectors at fixed time steps with 10 s intervals. This orbit file was generated per day and sent within 20 days of data acquisition. It was important to unveil the information on the data type, acquisition date, position, and orbit of the images during data acquisition.
These two images were selected to be combined in the image coregistration process to unify two unrelated data into a pair of data, i.e., master and slave images. The match between the two images could be observed from coherence value once the interferogram is formed. The entire pixel size of the slave image would be matched the pixel size of the master image. The image coregistration process would read two data. Each datum would be separated by its orbit to be replaced with an absolute orbit, i.e., the orbit of the satellite position when the satellite senses the surface of the earth.
The images with the correct orbital information were then transferred to the Back Geocoding process. Image resampling would be carried out using the Digital Elevation Model (DEM) which transformed the DEM coordinate system into an Earth-Centered Rotating (ECR) based Cartesian coordinate system. The resampling method was based on DEM and radar image resolution. Since the Sentinel-1 image has medium resolution and the DEM was derived from SRTM with WGS84 Geoid, it supported bilinear interpolation as a resampling method. The process of applying the satellite orbit and Back Geocoding which had geometric orientation was still insufficient for Sentinel-1 interferometric processing [4,6]. In order to increase the accuracy, Enhanced Spectral Diversity (ESD) was applied to estimate constant azimuth shifts between radar images by minimizing the phase discontinuity throughout the burst.
Interferometry was formed by multiplying the amplitude of the master image and the slave image when the phase difference of each image was displayed. The interferometric phase was generated from SAR image pixels which depended only on the difference between the flight path and the cell resolution in interferometry. In interferometric processing, elevation or deformation was estimated by eliminating other sources of errors. The flat-earth phase had to be removed from interferogram because it was produced by the curvature of the surface reference. This process could be done through the subtracted flat-earth phase, so that the topographic variations did not appear [3]. It was also intended to obtain accurate topographic or deformation information. The flat-earth phase was estimated using orbital information from metadata and subtracted from interferogram. The results obtained had coherent values (0.0 to 1.0/no useful information to perfect interferogram without noise). The minimum coherence value given by the European Space Agency (ESA) for the formation of the Digital Surface Height Model was 0.20. In order to be able to continue the process of estimating deformation using DInSAR required an image coregistration with a coherence value greater than 0.20. In this study, we mask out the coherence ≤ 0.5 in order to have a good quality interferogram. Low coherence presented in the interferograms could be caused by complex speckle noise from dense vegetation cover [46][47][48] and high rate of deformation [49][50][51].
The interferogram was produced from unified burst arrays. It needed an improvement in the form of merging all bursts into a complete interferogram through the TOPSAR Deburst process to identify range and azimuth. Each line in the direction of the sub-swath with the same time label would be merged with the adjacent sub-swaths. The interferogram could then be flattened by removing the topographic phase. DEM SRTM was used as a topographic model to perform the Topographic Phase Removal process. This step was implemented for deformation phase computation because eliminating topographic phase would remove the elevation value from the image pair interferogram and would leave only the deformation phase. The DEM SRTM acted only as topographic phase simulator, then it was removed from the interferogram of the processed image pairs.
In order to reduce noise in the fringes and to produce good deformation, the result from the previous step had to be filtered using Goldstein Phase Filtering by maintaining the fringes at the edges. The first process that occurs was an adaptive filter that was carried out before the frequency estimation was used to increase the accuracy of the estimation. Furthermore, to maintain the fringes characteristic, the estimated fringes frequency in each filtering patch was removed from the original noise phase [5,6]. Then the residual phase is refined using a Goldstein filter.
The DInSAR process produced images in radian unit (angular phase units) in the range −π to π, causing ambiguity problems. Although the deformation pattern could be observed, the main information regarding the deformation value could not be read properly. In order to get a deformation image containing metric value, an unwrapping process had to be carried out, and the absolute phase angle unit had to be changed into metric units in LOS (Line of Sight).
Unwrapping was one of the important steps to obtain high-quality phase unwrap results. The filtering process was done first, and the signal would be smoothed by exporting the subset results into the SANPHU format. Then the unwrapping phase was performed using the SNAPHU Unwrapping plugin. Image unwrapping was completed using the MCF (Minimum Cost Flow) method with the consideration that the processed image has a large area with a small coherence value [52][53][54][55].
A phase to displacement process was conducted to transform phase into deformation [4]. In order to find out the magnitude of deformation, a calculation of the displacement of the earth's surface was performed according to the equation along the line of sight (LOS) sensor [56,57]: where the value of λ is the wavelength of Sentinel-1 image, ∆∅ defo is the phase difference value and R is the wavelength distance. The result of deformation at this stage had not been georeferenced. Furthermore, the georeferencing process needed to be carried out so that the image had the correct coordinates on the map. Geocoding was performed by using DEM as a georeferenced model. In order to match the medium resolution of Sentinel-1, we applied bilinear interpolation for SRTM 1Sec HGT. Image coordinates were converted into the Geographic Coordinate System (WGS84). Indeed, radar coordinates were Range Doppler coordinates with some distortions, such as foreshortening and shadows, due to side-looking geometry. These coordinates should be converted into Geographic Coordinates System to match the coordinates of the Earth's surface. This process relied on Range Doppler Terrain Correction, in which we employed Range Doppler orthorectification method based on orbit data, radar recording time, slant range to the Earth's surface and referenced DEM. This produced final georeferenced deformation maps which contained uplift and subsidence value in meter. Methodological limitation of this study was caused by lack of time series terrestrial measurement. We were unable to acquire ground GPS points prior to each of the earthquake series (29 July; 5 August; 9 August; 19 August 2018), although our GPS survey after the disaster was successful. Therefore, we relied only on a satellite-based approach to achieve the aims of this study.

Magnitude of Displacement in Northern Part of Lombok Island
This study focused on earthquakes occurring on four dates: 29 July, 5 August, 9 August, and 19 August 2018 on northern Lombok Island. The spatial pattern of the ground displacement due to these earthquakes was visualized in an interferogram with radian units −π to π ( Figure 4) and a displacement map with metric values ( Figure 5). The interferogram of the foreshock earthquake of 6.4 Mw on 29 July 2018 was produced based on radar imageries on 27 July 2018 (master) and on 2 August 2018 (slave). Based on interferogram, the observed ground movement was focused on the northern shoreline. This involved a very light uplift of less than 0.155 m. The uplift was also found on the north-eastern shore (0.074 m). The subsidence was concentrated on the Rinjani crater rim, reaching −0.144 m.
The interferogram of the mainshock earthquake of 7.0 Mw on 5 August 2018 was derived from radar imageries on 2 August 2018 (master) and on 8 August 2018 (slave). Clear uplift was found from the north-western (0.574 m) to the north shore (0.452 m). Rinjani caldera was also uplifted less than 0.186 m. The uplift in the eastern shore was measured at 0.072 m. Light subsidence was observed at the western part of the island, near Mataram city. The interferogram of the aftershock earthquake of 6.2 Mw on 9 August 2018 was processed based on radar imageries on 8 August 2018 (master) and 14 August 2018 (slave). The subsidence was observed throughout the study area. It reached -0.217 m on the western shore and -0.389 m at the three touristic islands in front of the western shore (Gili Trawangan, Gili Meno and Gili Air island). Another subsidence of -0.130 m was also found at the Rinjani volcano.  We estimated displacements on a local scale: Mataram City, Pamenang, Tampes, Sukadana, Sembalun, and Belanting ( Figure 6). Mataram City appeared to be a modest average in subsidence and uplift, ranging from -0.053 m to 0.025 m. The most extreme displacements were found at Tampes, Pamenang and Sukadana. In Tampes, the highest uplift of 0.530 m was caused by a mainshock earthquake series of 7.0 Mw on 5 August 2018, while the deepest subsidence of -0.121 m was caused by the aftershock earthquake series of 6.2 Mw on 9 August 2018 ( Figure 6). Figure 5. Spatial distribution of vertical displacement due to earthquakes on 29 July, 5 August, 9 August, and 19 August 2018. Figure 5. Spatial distribution of vertical displacement due to earthquakes on 29 July, 5 August, 9 August, and 19 August 2018.

Average Magnitude of Displacement at Local Scale
We estimated displacements on a local scale: Mataram City, Pamenang, Tampes, Sukadana, Sembalun, and Belanting ( Figure 6). Mataram City appeared to be a modest average in subsidence and uplift, ranging from −0.053 m to 0.025 m. The most extreme displacements were found at Tampes, Pamenang and Sukadana. In Tampes, the highest uplift of 0.530 m was caused by a mainshock earthquake series of 7.0 Mw on 5 August 2018, while the deepest subsidence of −0.121 m was caused by the aftershock earthquake series of 6.2 Mw on 9 August 2018 ( Figure 6).
Pamenang was also affected by deformation since it experienced the deepest subsidence of −0.225 m due to the aftershock earthquake of 6. Most of the displacements that can be measured with good coherency are located next to the shoreline. Therefore, we conduct observation at 6 selected areas: Mataram City, Pamenang, Tampes, Sukadana, Sembalun, and Belanting. The highest uplift reaches 0.713 m at Tampes due to the mainshock earthquake of 7.0 Mw on 5 August 2018, while the deepest subsidence is measured -0.338 m at Pamenang following the aftershock earthquake of 6.2 Mw on 9 August 2018 ( Table 2). Mataram city has the modest displacement since it is located far from epicentrum.
Rapid assessment after event conducted by Pusat Gempa Nasional (Indonesian National Center for Earthquake studies) has revealed microatoll coral at northeastern part of the island [1] corresponding to an average uplift of 0.265 m in Belanting due to an earthquake of 6.9 Mw on 19 August 2018 that is measured in this study. Liquefaction is s also found at the subsidence area in Pamenang following the aftershock earthquake of 6.2 Mw on 9 August 2018 [1]. In the field, the deformation could also be corresponded to 4823 landslides after the earthquake on 5 August 2018 and 9319 landslides after the earthquake on 19 August 2018 [59] as well as small scale tsunamis following this earthquake series [60]. Dominant deformation at the northern area of Lombok Island indicates movement of Back Arc Trust in the north of the island [61].
The pattern and value of deformation in this research are similar to previous DInSAR research in the Lombok area, especially for earthquake on 5 August 2018 [52,4] and earthquake on 19 August 2018 [1,62]. However, earthquake on 29 July 2018 between this research and previous study [1] has shown different results because of the different master and slave images. This study applies radar images on 27 July 2018 (master) and on 2 August 2018, while PuSGen uses radar images on 25 July 2018 (master) and on 31 July 2018 (slave). The aftershock earthquake of 6.2 Mw on 9 August 2018 has never been analyzed before. Table 2. Zonal statistics of displacement due to earthquakes in selected areas (Mataram city, Pamenang, Tampes, Sukadana, Sembalun, and Belanting). Sembalun is located next to the Rinjani volcanic complex. It is a relatively flat area surrounded by mountainous area. Two major subsidences occurred due to the foreshock earthquake of 6.4 Mw on 29 July 2018 and the aftershock earthquake of 6.2 Mw on 9 August 2018. The subsidences were measured consecutively −0.095 m and −0.090 m. This was the only area having two nearly identical subsidence due to two different earthquakes.
Belanting showed an unique characteristic. The highest uplift due to the earthquake of 6.9 Mw on 19 August 2018 was found in this area. It reached 0.265 m. The deepest subsidence there was measured −0.110 m.

Discussion
Lombok earthquake series in 2018 contains four main events (29 July, 5 August, 9 August, and 19 August 2018) and at least 818 minor earthquakes in between these dates until 31 August 2018. This study emphasizes the displacement following each principal earthquake at the northern part of Lombok Island based on the DInSAR technique. The interferograms allow us to describe the pattern of displacement. The spatial distribution of vertical displacement in raster format helps to measure the value of uplift and subsidence in meters.
Each earthquake has different impacts on the displacement in Lombok Island. The first 6.4 Mw foreshock hitting on 29 July 2018 produces a minimum uplift effect on the island. However, it causes deep subsidences in the northwest area and on the Rinjani volcanic complex. The 7.0 Mw mainshock on 5 August 2018 causes an extreme uplift at the northern shore. The 6.2 Mw Aftershock on 9 August 2018 generates subsidence throughout the study area. The final earthquake of 6.9Mw on 19 August 2018 initiates a massive uplift in the study area and an extreme uplift at the northeastern shore. These could be caused by two reasons. Firstly, limited distance between Back Arc Thrust and the Lombok island produces greater deformation impacts at the northern shore ( Figure 2). Secondly, the epicenters are continuously changing between these four major earthquakes indicating unstable tectonic setting beneath the island due to brittle damage in the Earth's crust. Fortunately, there is no report on increasing volcanic activity. This confirms that there is no reciprocal relationship between volcanism and earthquake series [58].
Most of the displacements that can be measured with good coherency are located next to the shoreline. Therefore, we conduct observation at 6 selected areas: Mataram City, Pamenang, Tampes, Sukadana, Sembalun, and Belanting. The highest uplift reaches 0.713 m at Tampes due to the mainshock earthquake of 7.0 Mw on 5 August 2018, while the deepest subsidence is measured −0.338 m at Pamenang following the aftershock earthquake of 6.2 Mw on 9 August 2018 (Table 2). Mataram city has the modest displacement since it is located far from epicentrum. Rapid assessment after event conducted by Pusat Gempa Nasional (Indonesian National Center for Earthquake studies) has revealed microatoll coral at northeastern part of the island [1] corresponding to an average uplift of 0.265 m in Belanting due to an earthquake of 6.9 Mw on 19 August 2018 that is measured in this study. Liquefaction is s also found at the subsidence area in Pamenang following the aftershock earthquake of 6.2 Mw on 9 August 2018 [1]. In the field, the deformation could also be corresponded to 4823 landslides after the earthquake on 5 August 2018 and 9319 landslides after the earthquake on 19 August 2018 [59] as well as small scale tsunamis following this earthquake series [60]. Dominant deformation at the northern area of Lombok Island indicates movement of Back Arc Trust in the north of the island [61].
The pattern and value of deformation in this research are similar to previous DInSAR research in the Lombok area, especially for earthquake on 5 August 2018 [4,52] and earthquake on 19 August 2018 [1,62]. However, earthquake on 29 July 2018 between this research and previous study [1] has shown different results because of the different master and slave images. This study applies radar images on 27 July 2018 (master) and on 2 August 2018, while PuSGen uses radar images on 25 July 2018 (master) and on 31 July 2018 (slave). The aftershock earthquake of 6.2 Mw on 9 August 2018 has never been analyzed before.
Sentinel 1 imagery has been extensively used in global scale to analyze ground deformation due to single earthquakes [22][23][24][25][26]. DinSAR analysis derived from Sentinel 1 imagery on the Lombok earthquake series in Indonesia adds new insights. This small volcanic island with ultraplinian deposits and an active tectonic setting that was hit by multiple earthquakes has resulted in different ground deformation impacts on the spatio-temporal scale. The spatial distribution of uplift and subsidence has continuously changed over time because each earthquake produces an unique impact on ground deformation.

Conclusions
Lombok Island, Indonesia is located between the Back Arc Trust in the north and an active subduction zone of the Java trench. This tectonic setting leads to continuous exposure of this island to earthquakes. Recent shallow and great magnitude of earthquake series occurred in 2018. The first 6.4 Mw foreshock hit on 29 July 2018, followed by 7.0 Mw mainshock on 5 August 2018, then 6.2 Mw Aftershock on 9 August 2018, and finished by an earthquake of 6.9 Mw on 19 August 2018. This study discusses the displacement following each principal earthquake at the northern part of Lombok Island based on the DInSAR technique. The interferograms allow to describe the pattern of displacement. The spatial distribution of vertical displacement in raster format helps to measure the value of uplift and subsidence in meters. The highest uplift reaches 0.713 m at northern shore due to the 7.0 Mw earthquake on 5 August 2018, while the deepest subsidence is measured −0.338 m at the northwestern shore following the earthquake of 6.2 Mw on 9 August 2018. Dominant deformation at the northern area of Lombok Island indicates movement of Back Arc Trust in the north of the island. Insight from this study would be important to local authorities to evaluate existing earthquake's impacts and prepare the mitigation strategies to face earthquake induced ground displacement.