Identification of Sea Level Rise and Land Subsidence Based on Sentinel 1 Data in the Coastal City of Pekalongan, Central Java, Indonesia

Sea level rise is a pure impact of climate change. However, the process of studying sea level rise must include local factors that influence such as land subsidence. This study focuses on sea level rise using the CEEMDAN method and land subsidence using the DInSAR method. The location of this research is Pekalongan, Central Java, Indonesia. Tidal data used in this study was for five years, from 2016 to 2020, obtained from the Geospatial Information Agency (BIG). Then the data used to study land subsidence in this study uses Sentinel-1 Synthetic Aperture Radar (SAR) data in 2015, 2016, 2017, and 2020. Pekalongan is an area with mixed diurnal tidal types with Formzahl number 1.7. The sea level rise in Pekalongan is relatively high, at 10.6 mm/year. Then the land subsidence that occurred in Pekalongan is the phenomenon that has the most influence on the occurrence of coastal flooding in the region. The average land subsidence on the coast of Pekalongan is 5.37 cm/year. In addition, the sampling results in 6 areas showed that the most significant decrease was in Area 2, with a decrease of -7.91 cm/year. Based on this research, land subsidence is the most considerable influence on flooding in Pekalongan compared to sea level rise.


INTRODUCTION
Climate change that occurs on a global scale impacts the lives and activities of coastal communities. One threat that occurs due to climate change to coastal areas is the damage to coastal ecosystems due to the inundation of coastal areas due to sea level rise (SLR) (Mehyar et al., 2018). Sea level rise is generally understood as a pure impact of climate change. However, studying sea level rise must include local factors that influence this phenomenon (Suroso and Firman, 2018). Based on USACE (2014), one climate change factor alone will not significantly impact. Therefore, it is necessary to study several factors of climate change to determine the very severe impact of climate change.
Sea level rise is divided into relative sea level rise (relative SLR) and global geocentric sea level rise (global geocentric SLR) (Bott et al., 2021). Relative SLR is a sea level rise condition observed directly along the coastline. Relative SLR is strongly influenced by land subsidence, which is a significant threat in coastal areas (Esteban et al., 2020). Meanwhile, global geocentric SLR is an SLR observation from an altimetry radar satellite. So, the results of the SLR shown do not consider the effect of land subsidence (Mehyar et al., 2018).
The northern region of Java Island is a coastal area that is very vulnerable to sea level rise, especially those influenced by land subsidence (Sarah et al., 2021). Nerem et al. (2018) and Esteban et al. (2020) have researched the study of sea level rise, relative and global geocentric, in the northern region of Java Island, precisely north of Jakarta and Semarang. The results show that relative SLR due to land subsidence in northern Jakarta has a value of 12 cm/year, while global geocentric SLR has a value of -3.2 mm/year. Relative SLR in the north of Semarang shows a 10 cm/year value, which is based on the GNSS-corrected tide gauge in Semarang (Niesters et al., 2018). The previous research shows that the relative value of SLR due to land subsidence in the north of Java Island is greater and significantly impacts coastal areas. Thus, it is essential to study the factors of land subsidence on relative sea level rise (relative SLR) on the north coast of Java Island.
Besides Jakarta and Semarang, one of the cities on the north coast of Java Island, which is also affected by sea level rise due to land subsidence, is Pekalongan City. Andreas et al. (2018) has researched land subsidence in Pekalongan City, which shows a very high value of 10-14 cm/year. This value is equivalent to the highest level of land subsidence in Jakarta and Semarang. In addition, (Adi Iskandar et al., nd) conducted a sea level rise study with one month of tidal data in February 2020 using a conventional method in the form of linear regression showing that the sea level rise that occurred in Pekalongan was around 4.3 mm/year. Studies on sea level rise by considering land subsidence have been carried out. Remote sensing and GIS analysis are the most widely used methods to assess land subsidence. Yu et al. (2022) have conducted a spatial analysis in coastal areas using Synthetic Aperture Radar (SAR) for coastal dynamics, which shows changes in coastline due to land subsidence. Zainuri et al. (2022) also researched the land subsidence on relative sea level rise in Pekalongan City using Sentinel-1 Synthetic Aperture Radar (SAR) with Single Band Algorithm (SBA) differential interferometry and satellite altimetry tide data. The study shows an increasing trend in tidal flood inundation areas due to land subsidence and sea level rise every year. Nashrrullah et al. (2013) conducted a land subsidence study using ALOS/PALSAR data from 2008 to 2009 and found that the land subsidence in Pekalongan was 3 cm/year.
In addition to level land subsidence analysis, tidal data analysis to assess and identify sea level rise is also carried out. Tidal data analysis can be performed using spectral analysis to decompose non-linear and non-stationary signals to obtain rising sea-level trends. Non-linear analysis of tidal trends can be performed using the Complete Ensemble Empirical Mode Decomposition (CEEMD) and Hilbert-Huang transformation (HHT) methods. Anwar et al. (2022) have conducted a study to analyze tidal trends using this method. Sea level elevation data is decomposed using HHT based on empirical mode decomposition (EMD). Breaker and Ruzmaikin (2013) has also researched to analyze the acceleration of sea level rise by applying ensemble empirical mode decomposition (EEMD). In EEMD, the analysis is done by adding noise in the decomposition to stabilize the process (Han and van der Baan, 2013). Previous studies have studied sea level rise influenced by land subsidence. However, applying the spectral analysis method in processing tidal data in the Pekalongan City area has yet to be carried out. This study aims to identify sea level rise by considering land subsidence. Spatial analysis was conducted to study the subsidence in the study area using Synthetic Aperture Radar data with the DInSAR (Differential Interferometric Synthetic Aperture Radar) method. The trend of sea level rise was analyzed using the complete ensemble empirical mode decomposition with adaptive noise (CEEMDAN) spectral analysis method.

MATERIAL AND METHOD
In order to study sea level rise in the Pekalongan region, Central Java, Indonesia (Figure 1), this study uses tidal data recorded by the Geospatial Information Agency (BIG) via fixed stations. The tidal data used is five years, from 2016 to 2020. The period of the tidal dataset used in this study is 1 hour (Figure 2).
Then to perform land subsidence analysis, this study used Synthetic Aperture Radar (SAR) data from Sentinel-1. Available data are from 2015, 2016, 2017, and 2020. SAR is a system capable of obtaining complex, high-resolution images of extensive areas of terrain, usually located on orbital or airborne platforms. However, it can also be used in ground-based deployments (Yerro et al., 2014). Detailed characteristics of the downloaded data are listed in the following list.
Sentinel-1 image data was processed using the open-source Sentinel Applications Platform (SNAP) software. The open-source policy and excellent software capabilities make SNAP the main choice for processing SAR data and optical image processing (Tzouvaras et al., 2019). Then, the SNAPHU software was used for the unwrapping process of the methodology.

Tidal Analysis
Data were analyzed for harmonics using T_TIDE to obtain 68 tidal harmonic components. These results can be used to analyze sea level data and determine the residual sea level that is influenced other than tides, such as the result of meteorological influences. Pawlowicz et al. (2002) used an algorithm to estimate the amplitude and phase based on the FORTRAN algorithm developed by Godin (1972), Foreman (1977), and Foreman (1978). However, in contrast to the FORTRAN algorithm, Pawlowicz et al. (2002) uses complex algebra directly instead of using the sine and cosine fit separately. Then, identify the type of tide using form factor F as follows.

= ( 1 + 1)/( 2 + 2)
For F<0.25 semidiurnal tides, F between 0.25 and 1.5 mixed semidiurnal tides, F between 1.5 and 3.0 mixed diurnal tides, and F is more than 3.0 diurnal tides. Wu et al. (2007) define the trend function as an intrinsically fitted monotonic function or function with at most one extreme value in a specific data range. The Empirical Mode Decomposition (EMD) method adopted can maintain the fundamental nature of nonlinear and non-stationary time series data, that the appropriate physical interpretation of the trend function in each data range should not change with the addition of new data (Wu et al., 2007). EMD decomposes the time series data into a finite set of signals called Intrinsic Mode Functions (IMFs). IMFs represent different oscillations in the data (Han and van der Baan, 2013). The EMD method is non-parametric. The trend function can represent the low-frequency variability of the sea level whose period is longer than the data range, and this will be useful for investigating how the sea level trend varies (Chen et al., 2014).

Complete Ensemble Empirical Mode Decomposition with Adaptive Noise (CEEMDAN)
EMD has shortcomings in its mixing mode. Therefore, Wu and Huang (2007) proposed a new method, namely Ensemble Empirical Mode Decomposition, as an adaptive time-frequency data analysis method. EEMD can decompose the noise ensemble copy of the original signal. So that the result can be obtained using an average. The addition of white Gaussian noise can reduce the mixing mode to some extent. However, according to Gang et al. (2017), this causes new problems, the difference in signal realizations plus noise that allows for discrepancies in the mode of the number, thus complicating final averaging. Therefore, Colominas et al. (2014) proposed the Complete Ensemble Empirical Mode Decomposition with Adaptive Noise (CEEMDAN) algorithm, which has proven to be a significant development of EEMD. This method can ignore reconstruction errors and solve the problem of the number of modes that are different for different signals plus noise realization. The CEEMDAN algorithm can be explained as follows.
Step 1: Based on the main signal x(t), xi(t) is obtained using the following equation.
( ) = ( ) + 0 [ ( )] = 1,2,3 … where ( ) shows the realization of white Gauss noise with zero mean and unit variance, is the operator that yields the main mode obtained from EMD, 0 is a constant.
Step 2: Identify all the extrema points (local maxima and minima) of the ( ) and interpolate them (cubic spline interpolation) to find the upper and lower envelopes. Then calculate the mean envelope [ ( )] to find the first residual 1 ( ): Step 3: Calculate the first mode with the following equation. Step 5: For = 3, 4 … repeat Step 2 -Step 4, and can find the residual k th : Step 6: Calculating the k th : Step 7: Finally, the original signal x(t) can be reconstructed using the following equation.

Land subsidence
This study uses the DInSAR method to determine how much the land surface changes in the coastal area of Pekalongan. The DInSAR method is used to see changes in the baseline in SAR data so that the phase difference can later be reduced to a change, meaning that the information contained in the interferometric phase is calculated as the phase difference between two SAR images obtained at two different times but at almost the exact location (Fajrin et al., 2021). DInSAR can be processed using SNAP software developed by the European Space Agency (ESA) (Mandal et al., 2019). When the interferometric phase appears as follows: Note: ∆ = the contribution of the flat earth phase; ∆ = the topography; ℎ ℎ = the land deformation measured along the Line of Sight (LOS); the ℎ = the contribution of the atmosphere to the interferometric phase and is generic noise.
The process starts with data collection, and then it is estimated to see the relationship between the master and slave image pairs between 2015 and 2016, 2016 and 2017, and 2017 and 2020. After getting the master and slave pair, the slave image is adjusted to the master image during the image correction section. Based on the co-registered image, this section forms an interferogram and estimates the coherence of the master and slave image pairs. The next stage is to make interferogram and coherence. The lower the consistency ( < 0.2), the lower the level of adjustment between these images (Bamler and Hartl, 1998). If the coherence is greater than 0.2, the step can be continued to the TOPSAR Deburst stage. At this stage, the black line is separated between the bursts. The image cropping section to focus on this research is carried out after the TOPSAR Deburst process. The next step is to eliminate the topographic phase, where in this section, create an interferogram simulation based on the Digital Elevation Model (DEM) and divide the topographic phase of the process interferogram. The DEM used in this study is the SRTM 3-sec DEM which is automatically downloaded to the SNAP software.
Creating an interferogram requires filtering to reduce phase noise. The filter process uses the Goldstein Filter, available in the SNAP feature. After filtering, the next step is Multilooking with a function to reduce phase noise in SAR images by forming pixels resembling squares. However, the phase formed in the Multilooking process still has an ambiguous phase, so removing the ambiguous phase into an absolute phase is necessary. This stage is called the Unwrapping Phase, which is processed using the SNAPHU software. The next stage is to convert the unwrap data into phase to height to determine the height difference from the DInSAR process or convert from slant to height using Phase to Displacement in SNAP software. The next step is Range Doppler Terrain Correction to put the image in its actual position. The processing results are then exported into NetCDF form (.nc) for layout using Python.

RESULT AND DISCUSSION
Based on the results of tidal processing using T_TIDE, obtained 68 tidal harmonic components as follows (Table 2). Based on the calculation results, the Formzahl number is 1.82, which indicates that the area has a mixed diurnal tidal type. This type can be confirmed by the tides' pattern forming a mixed diurnal. In addition, Sedyoko et al. (2013) conducted a tidal study in the Pekalongan area with the results of the Formzahl number 1.7, which indicates a mixed diurnal tidal type. In this study, we apply CEEMDAN to tidal data for five years (2016)(2017)(2018)(2019)(2020). In this case, the effect of annual and semiannual cycles is not removed because it will separate from the long-term tidal trend on the non-lower decomposition method. Annual trends and semiannual cycles will be separated into separate IMF results. In this study, separate IMF 10 -IMF 12. Figure 2 shows that there are 12 IMFs and residuals consisting of mean sea level, steric sea level, the ocean mass, in descending order of the high-frequency noise, the annual cycle, the interannual variability, and the trend function. When CEEMDAN decomposes interannual variability with different time scales into separate IMFs, it combines IMFs with time scales longer than the annual cycle and shorter than the last trend function to produce IMFs of interannual variability. This research focuses on the latest IMF results in the form of a sea level rise trend shown on the y-axis (Figure 2). The y-axis shows units of centimeters. It will be converted to millimeters when calculating the sea level rise value.
The thirteenth graph (Figure 3) shows a trend function of the mean sea level tie in the Pekalongan area from 2016 to 2020. According to Wu et al. (2007), the trend function depends on the length of the time series data. This trend is the latest mode which shows an incomplete wave so that the CEEMDAN algorithm does not recompose it. Based on the residual graph, the trend significantly increased from 2016 to 2020. IMF 10 -IMF 12, an annual and semiannual cycle, shows that, in the end, it is affected by climatological conditions, namely ENSO (El Nino Southern Oscillation). At the IMF, there is an increase at the end of the time series data. However, for the mean sea level and mean steric sea level, the impact of this ENSO does not change the fundamental nature of the trend function, especially how ENSO varies during the study year (Chen et al., 2014). In Figure 4 shows a trend in the residuals, while Figure 3 on a larger scale. Figure 4 shows an upward trend from 146 cm to 151 cm during five years. So, the sea level rise due to climate change in the Pekalongan area is 10.6 mm/year. This increase is relatively significant compared to the Global Mean Sea Level studied by Chen et al. (2014), which experienced a sea level rise of 3.2 mm/year. Figure 5 shows the distribution of land subsidence in the Pekalongan region. It shows that the center of subsidence is on the west side of the coast of Pekalongan. Meanwhile, the east side shows a significant land increment. Based on the distribution map, the land subsidence occurred in central Pekalongan City. On the coastal side of Pekalongan, the average land subsidence is -5.37 cm/year.
In this study, the average land subsidence in the Pekalongan area was sampled in each area until it was divided into six areas, as shown in Figure 5. This is done to determine the variation of subsidence or increment around the coast of Pekalongan. Sampling is also carried out on annual land subsidence maps so that variations in subsidence or increment can be seen based on the year at the 6 points.  Based on Table 3, Area 2 is the area that experienced the most decline, with a very high decline. The average decline in Area 2 reached -7.91 cm/year. Then, the area with the lowest land subsidence value is Area 6, with -2.622 cm/year. Based on the results, the west coast of Pekalongan experienced more land subsidence than the east coast of Pekalongan.
The Pekalongan area is dominated by land subsidence, shown in Table 3, in which six areas experienced land subsidence. This aligns with research conducted by Ardha et al. (2022). The Pekalongan area experienced land subsidence of 11.2 cm/year. Pekalongan is located in a smooth alluvial area, and population development is significant; this causes land subsidence in the Pekalongan area to be quite large. Nashrrullah et al. (2013) also revealed that land subsidence in the Pekalongan area was 4.8 cm/year. Like other areas, the land in the Pekalongan area is alluvial soil that has yet to be consolidated and compacted perfectly. However, settlements on this land have been very rapid, coupled with the growing industry in Pekalongan, which has resulted in significant land subsidence. This significant subsidence has a major impact on the morphological conditions of the surrounding coast and has an impact on the occurrence of flooding in Pekalongan. So that the flooding in Pekalongan is influenced mainly by land subsidence because it has a huge number compared to sea level rise.

CONCLUSIONS
In this study, we studied the sea level rise using the latest spectral method, namely the Hilbert Huang Transform (HHT) with the Complete Ensemble Empirical Mode Decomposition with Adaptive Noise (CEEMDAN) decomposition method and conducted a study on land subsidence using the Differential Interferometric Synthetic Aperture Radar (DInSAR) method. Tidal data was used from 2016 to 2020, and landside data used are 2015, 2016, 2017, and 2020. Based on the research, the Pekalongan area's tidal type is mixed diurnal with a Formzahl value of 1.7. Then based on the results of the CEEMDAN decomposition, there are 12 Intrinsic Mode Functions (IMFs) and one residual. IMF 10-IMF 12 are annual and semiannual cycles that affect the tides. Then the residual is the sea level rise trend obtained from the decomposition results. The sea level rise in Pekalongan is 10.6 mm/year. This value is quite significant for the sea level rise category. Then for land subsidence based on the average sampling results in 6 coastal areas of Pekalongan, some areas have the most severe land subsidence, namely Area 2, with a decrease of -7.9 cm/year. Based on this, the effect of land subsidence is greater than the sea level rise, which results in flooding in the coastal area of Pekalongan. This study can be used as a reference for further researching the occurrence of sea level rise and land subsidence in the Pekalongan area by adding other supporting data such as altimetry satellite data.