Shoreline Dynamics in East Java Province, Indonesia, from 2000 to 2019 Using Multi-Sensor Remote Sensing Data

: Coastal regions are one of the most vulnerable areas to the effects of global warming, which is accompanied by an increase in mean sea level and changing shoreline conﬁgurations. In Indonesia, the socioeconomic importance of coastal regions where the most populated cities are located is high. However, shoreline changes in Indonesia are relatively understudied. In particular, detailed monitoring with remote sensing data is lacking despite the abundance of datasets and the availability of easily accessible cloud computing platforms such as the Google Earth Engine that are able to perform multi-temporal and multi-sensor mapping. Our study aimed to assess shoreline changes in East Java Province Indonesia from 2000 to 2019 using variables derived from a multi-sensor combination of optical remote sensing data (Landsat-7 ETM and Landsat-8 OLI) and radar data (ALOS Palsar and Sentinel-1 data). Random forest and GMO maximum entropy (GMO-Maxent) accuracy was assessed for the classiﬁcation of land and water, and the land polygons from the best algorithm were used for deriving shorelines. In addition, shoreline changes were quantiﬁed using Digital Shoreline Analysis System (DSAS). Our results showed that coastal accretion is more profound than coastal erosion in East Java Province with average rates of change of +4.12 (end point rate, EPR) and +4.26 m/year (weighted linear rate, WLR) from 2000 to 2019. In addition, some parts of the shorelines in the study area experienced massive changes, especially in the deltas of the Bengawan Solo and Brantas/Porong river with rates of change (EPR) between − 87.44 to +89.65 and − 18.98 to +111.75 m/year, respectively. In the study areas, coastal erosion happened mostly in the mangrove and aquaculture areas, while the accreted areas were used mostly as aquaculture and mangrove areas. The massive shoreline changes in this area require better monitoring to mitigate the potential risks of coastal erosion and to better manage coastal sedimentation.


Introduction
Global warming caused by increasing greenhouse gas abundance in the atmosphere is projected to result in an average sea level rise of 20 to 200 cm in the 21st century [1], and the sea level in the year 2100 is projected to be 310 ± 30 mm higher than that in 1990 [2]. The sea level rise induced by global warming is thought to be accelerated by the contribution of ice sheet melting in Greenland and Antarctica [3][4][5]. The increase in sea level has the potential to affect shorelines [6] and thus to amplify coastal erosion [7][8][9][10]. In addition, global warming increases the risk of extreme weather by strengthening the thermodynamic factors for generating tropical cyclone and storm surges, which subsequently increases the risk of coastal flooding and coastal erosion due to the elevated surge height from the intensifying winds, waves, and atmospheric pressures [11][12][13][14]. These climate-changeinduced events could potentially affect communities living near coasts, meaning that global shoreline changes in this area are also relatively understudied because of a lack of big datasets necessary for regional-level mapping. Therefore, our study aimed to 1. Implement the cloud computing platform of the Google Earth Engine and its accompanying remote sensing datasets (i.e., Landsat-5 TM, Landsat-7 ETM, Landsat-8 OLI, and Sentinel-1) to perform a regional mapping of the shoreline dynamics in a part of East Java Province; 2.
Monitor the regional shoreline dynamics and the affected land use in a part of the eastern coastal areas of East Java Province from 2000 to 2019 at 4-to 5-year intervals.

Study Area
The study took place in the eastern coastal region of East Java Province, Indonesia, between 112.45 to 114.21º E and 6.833 to 7.93º S (Figure 1). East Java Province was chosen as the study area because of the presence of Surabaya, a coastal megacity, and the industrial development of the coastal regions. This is compounded by the degradation of the watersheds of the main rivers in East Java (i.e., Brantas/Porong and Bengawan Solo), which has increased organic sediment load [39] and domestic plastic pollution [40,41] due to improper waste management, both of which affect the marine environment [42].
Land 2020, 9, x FOR PEER REVIEW 3 of 16 shoreline changes in this area are also relatively understudied because of a lack of big datasets necessary for regional-level mapping. Therefore, our study aimed to 1. Implement the cloud computing platform of the Google Earth Engine and its accompanying remote sensing datasets (i.e., Landsat-5 TM, Landsat-7 ETM, Landsat-8 OLI, and Sentinel-1) to perform a regional mapping of the shoreline dynamics in a part of East Java Province; 2. Monitor the regional shoreline dynamics and the affected land use in a part of the eastern coastal areas of East Java Province from 2000 to 2019 at 4-to 5-year intervals.

Study Area
The study took place in the eastern coastal region of East Java Province, Indonesia, between 112.45 to 114.21ºE and 6.833 to 7.93ºS ( Figure 1). East Java Province was chosen as the study area because of the presence of Surabaya, a coastal megacity, and the industrial development of the coastal regions. This is compounded by the degradation of the watersheds of the main rivers in East Java (i.e., Brantas/Porong and Bengawan Solo), which has increased organic sediment load [39] and domestic plastic pollution [40,41] due to improper waste management, both of which affect the marine environment [42].

Data
This study used Landsat surface reflectance data with a 30 m spatial resolution as the main datasets, which were accessed and processed in the GEE for the years 2000, 2005, 2010, 2015, and 2019 by taking the median reflectance values within each year. Depending on availability, some years (i.e., 2010, 2015, and 2019) included SAR data from ALOS Palsar/Global Palsar (horizontal-horizontal, HH, and horizontal-vertical, HV, polarizations at a 25 m spatial resolution) and Sentinel-1 C-band data (vertical-vertical, VV, and vertical-horizontal, VH, polarizations at a 10 m spatial resolution) to complement the Landsat surface reflectance (SR) data in the GEE. All radar data were resampled to a 30 m spatial resolution to match the resolution of the Landsat data. The combination of SAR backscatter and spectral bands from optical data is thought to increase the accuracy of land-use

Data
This study used Landsat surface reflectance data with a 30 m spatial resolution as the main datasets, which were accessed and processed in the GEE for the years 2000, 2005, 2010, 2015, and 2019 by taking the median reflectance values within each year. Depending on availability, some years (i.e., 2010, 2015, and 2019) included SAR data from ALOS Palsar/Global Palsar (horizontal-horizontal, HH, and horizontal-vertical, HV, polarizations at a 25 m spatial resolution) and Sentinel-1 C-band data (vertical-vertical, VV, and vertical-horizontal, VH, polarizations at a 10 m spatial resolution) to complement the Landsat surface reflectance (SR) data in the GEE. All radar data were resampled to a 30 m spatial resolution to match the resolution of the Landsat data. The combination of SAR backscatter and spectral bands from optical data is thought to increase the accuracy of land-use and land-cover classification, as reported in studies on classifying land cover in Southeast Asia [43], identifying flooded wetland areas [44], and enhancing the accuracy of shoreline mapping [45].
Several spectral indices from the Landsat data were generated as additional variables, besides the spectral bands, for the input analysis. Several indices were used, such as the normalized difference water index (NDWI) [46] (Equation (1)), the normalized difference vegetation index (NDVI) [47] (Equation (2)), and the non-shadow version (AWEInsh) of the automated water extraction index (AWEI) [48] (Equation (3)). The non-shadow version was used because of the relatively clear shadow pixels (i.e., cloud or dark surfaces) in the study areas. The equations for these indices are provided below.
where blue, green, red, near infrared (NIR), shortwave infrared 1 (SWIR1), and shortwave infrared 2 (SWIR2) represent the spectra of the optical bands of the Landsat data. The complete list of variables used for classification can be found in Table 1.

Supervised Classification
The supervised classification in this study was performed using GMO maximum entropy (GMO-Maxent) in the GEE [49]. GMO-Maxent is a multinomial logistic regression model for multi-class classification that is less computationally intensive than the standard optimization of the conditional maximum entropy (Maxent) model for deriving the probability distribution. The GMO-Maxent algorithm is rooted in the concept of entropy or amount of information stored within a variable [50,51]. The entropy was calculated from the sum of the occurrence probabilities of the variables p(i) multiplied by the surprise element or the binary logarithm (log 2 p(i)) of the variable i. Therefore, Shannon's entropy (E) can be summarized as Equation (4): The concept of Shannon's entropy was brought into its modern form, usually known as the maximum entropy principle, by Jaynes [52]. This form was developed to deal with the unknown estimation of the probabilities of an event p(i) by using the probability distribution of the maximum entropy from a known event, as in statistical mechanics. The maximum entropy probability distribution can be viewed as the most uniform and spread-out distribution that can satisfy any given constraint to represent the target distribution [53,54].
This algorithm is rarely used for the classification of remote sensing data [55], especially for multi-class classification. However, it is a well-known method in studies of species distributions and habitat models [53,56,57], and it has been reported to perform better in one-class or binary classification applications than other algorithms, such as the support vector machine method [58]. The GMO-Maxent method was compared with the random forest (RF) method [59] in classifying land, suspended sediment, and seawater body classes. To assess the accuracy of both methods, we randomly distributed 500 points, both on land and non-land bodies, along 1 km of the detected shoreline and constructed a confusion matrix to calculate the producer's accuracy (PA), user's accuracy (UA), and overall accuracy (OA) ( Figure 2).
Land 2020, 9, x FOR PEER REVIEW maximum entropy probability distribution can be viewed as the most unifor spread-out distribution that can satisfy any given constraint to represent the target bution [53,54].
This algorithm is rarely used for the classification of remote sensing data [55 cially for multi-class classification. However, it is a well-known method in studies cies distributions and habitat models [53,56,57], and it has been reported to perform in one-class or binary classification applications than other algorithms, such as the s vector machine method [58]. The GMO-Maxent method was compared with the r forest (RF) method [59] in classifying land, suspended sediment, and seawater bod ses. To assess the accuracy of both methods, we randomly distributed 500 points, b land and non-land bodies, along 1 km of the detected shoreline and constructed a sion matrix to calculate the producer's accuracy (PA), user's accuracy (UA), and accuracy (OA) ( Figure 2). In addition to the land and seawater classification, additional supervised cla tion with RFs was used to identify the affected areas and the type of land use in 2019 from the results of coastal erosion and accretion. Two areas with major change identified from the shoreline analysis, and supervised classification was conducte following classes were identified from the imagery: (1) aquaculture, (2) mangrove, ricultural, (4) built-up, and (5) bare land areas.

Shoreline Refinement and Change Assessment
The land polygons of each classified map from 2000 to 2019 were used to rep the shorelines and to quantify shoreline changes. After selecting the best results fr classification and before quantifying the shoreline changes, visual and manual refin were conducted by identifying any misclassifications that overshot or undershot la ter boundaries and by editing the vector vertices to better represent the border areas. This refinement was conducted to ensure that the calculations would mos rately reflect the shoreline changes. In addition to the land and seawater classification, additional supervised classification with RFs was used to identify the affected areas and the type of land use in 2000 to 2019 from the results of coastal erosion and accretion. Two areas with major changes were identified from the shoreline analysis, and supervised classification was conducted. The following classes were identified from the imagery: (1) aquaculture, (2) mangrove, (3) agricultural, (4) built-up, and (5) bare land areas.

Shoreline Refinement and Change Assessment
The land polygons of each classified map from 2000 to 2019 were used to represent the shorelines and to quantify shoreline changes. After selecting the best results from the classification and before quantifying the shoreline changes, visual and manual refinement were conducted by identifying any misclassifications that overshot or undershot land/water boundaries and by editing the vector vertices to better represent the border of land areas. This refinement was conducted to ensure that the calculations would most accurately reflect the shoreline changes.
To estimate the rate of shoreline changes, the ArcGIS plugin of Digital Shoreline Analysis System (DSAS) version 5.0 [60] was used. This plugin has been employed in various studies to calculate shoreline changes. It provides various statistical analyses, including linear regression rate (LRR), end point rate (EPR), and weighted linear regression (WLR) [61,62], as well as net shoreline movement (NSM) and shoreline change envelope (SCE) estimation [60]. EPR provides the rate of change based on the distance difference between the shoreline position in the baseline year and the end year, while WLR provides the best-fitted line for the linear rate of change by considering the uncertainty value [32]. EPR can yield a negative value, which represents coastal erosion, and a positive value, which represents coastal accretion. DSAS provides the capability for detailed and systematic calculation of the shoreline change rate by automatically generating transect lines perpendicular to a baseline [10,60]. The statistical analysis in DSAS was carried out by placing the transect lines along the shoreline, with the starting points of the transect lines being determined by the baseline. The baseline that we used in this study was the shoreline in the year 2000. In addition, we defined the following parameters for DSAS: a transect spacing of 50 m and a maximum search distance of 2 km.
A summary of the analysis workflow is presented in Figure 3. The GEE code for conducting the analysis is provided in the Supplementary Materials. To estimate the rate of shoreline changes, the ArcGIS plugin of Digital Shoreline Analysis System (DSAS) version 5.0 [60] was used. This plugin has been employed in various studies to calculate shoreline changes. It provides various statistical analyses, including linear regression rate (LRR), end point rate (EPR), and weighted linear regression (WLR) [61,62], as well as net shoreline movement (NSM) and shoreline change envelope (SCE) estimation [60]. EPR provides the rate of change based on the distance difference between the shoreline position in the baseline year and the end year, while WLR provides the best-fitted line for the linear rate of change by considering the uncertainty value [32]. EPR can yield a negative value, which represents coastal erosion, and a positive value, which represents coastal accretion. DSAS provides the capability for detailed and systematic calculation of the shoreline change rate by automatically generating transect lines perpendicular to a baseline [10,60]. The statistical analysis in DSAS was carried out by placing the transect lines along the shoreline, with the starting points of the transect lines being determined by the baseline. The baseline that we used in this study was the shoreline in the year 2000. In addition, we defined the following parameters for DSAS: a transect spacing of 50 m and a maximum search distance of 2 km.
A summary of the analysis workflow is presented in Figure 3. The GEE code for conducting the analysis is provided in the Supplementary Materials.

Accuracy Assessment
The accuracy assessment of RFs and GMO-Maxent showed that GMO-Maxent performed best in separating land and water body classes. Table 2

Accuracy Assessment
The accuracy assessment of RFs and GMO-Maxent showed that GMO-Maxent performed best in separating land and water body classes. Table 2  The better accuracy of RFs in 2015 and 2019 can be attributed to the use of Landsat-8 and Sentinel-1 data in the classification inputs, which helped to increase the accuracy of the RF algorithm. However, the other multi-sensor combination of Landsat-7 ETM and Global PALSAR-2 yearly mosaic data in 2010 did not yield the same results, with the accuracy of the RF models still falling below that of the GMO-Maxent models. It should be noted that Landsat-7 ETM in 2010 suffered from a scan-line error, rendering missing values on several lines in the imagery, which were filled by annual median values, but this was still imperfect, with noticeable changes in the spectral values of the surrounding pixels. It seems that the RF method was affected more by this artifact than the GMO-Maxent model. For the years 2000 and 2005, the RF algorithms yielded 6% to 9.6% lower overall accuracies than GMO-Maxent. The low accuracies of the RF method for these years resulted from the large number of false positives (commission errors), as evidenced by the fact that the user's accuracies were lower than the producer's accuracies in the land class. These errors can be attributed largely to the presence of fishponds in some parts of the study area, which were classified as bodies of water and thus affected the shoreline results. Therefore, if the RF results were used, it would require laborious refinements to be able to accurately estimate the shoreline changes for those areas (Figure 4a,b). In conclusion, the shoreline change analysis was ultimately conducted using the results of the GMO-Maxent method, which yielded a more stable classification with relatively few errors, thus reducing the need for manual refinement.

Shoreline Changes from 2000 to 2019
The Although the statistics of average EPR and WLR indicated overall coastal accretion, we identified regions of both major erosion and major accretion ( Figure 6A) in the northwestern part of East Java Province, which represents the delta of the Bengawan Solo river (Figure 6A), and coastal accretion, which represents the delta of the Brantas/Porong river (Figure 6B).  Although the statistics of average EPR and WLR indicated overall coastal accretion, we identified regions of both major erosion and major accretion ( Figure 6A) in the northwest-ern part of East Java Province, which represents the delta of the Bengawan Solo river ( Figure 6A), and coastal accretion, which represents the delta of the Brantas/Porong river ( Figure 6B).   (Figure 7). These results, while indicating the overall dominance of sedimentation in the area, also point to the delta's highly dynamic nature. The land covering this delta is dominated by aquaculture ponds, in ad-   For the Brantas/Porong delta region ( Figure 6B), we found an average shoreline change of +13.28 m/year from 2000 to 2019. The rate of shoreline change ranged from −19.98 to +111.75 m/year (EPR) (Figure 8). Unlike the Bengawan Solo delta, the Brantas/Porong delta showed no erosion at rates of more than ~20 m/year. The land-cover composition in the Brantas/Porong delta is similar to that in the Bengawan Solo delta; the acquired land is used for aquaculture ponds and as mangrove conservation areas. The massive coastal sedimentation in the delta of the Brantas/Porong river is due to a mud volcano explosion in May 2006. During this event, underground mud was lifted to the surface, carried by the river, and deposited in the delta, adding around 3,881.8 ha of additional land to the area [63].

Land Use of the Accretion and Erosion Areas in the Coastal Regions in 2000 to 2019
The mapped land-use classes in the Bengawan Solo and Brantas/Porong deltas can be found in Figure 9. Most of the delta regions were used as mangrove and aquaculture areas, especially in the Brantas/Porong delta. In contrast, in the Bengawan Solo delta, small portions of the areas were used as plantation and built-up areas. Based on the landuse classification results and the shoreline changes in 2000 and 2019, the mangrove and aquaculture areas were the main affected land-use classes in those deltas.  (Figure 8). Unlike the Bengawan Solo delta, the Brantas/Porong delta showed no erosion at rates of more than~20 m/year. The land-cover composition in the Brantas/Porong delta is similar to that in the Bengawan Solo delta; the acquired land is used for aquaculture ponds and as mangrove conservation areas. The massive coastal sedimentation in the delta of the Brantas/Porong river is due to a mud volcano explosion in May 2006. During this event, underground mud was lifted to the surface, carried by the river, and deposited in the delta, adding around 3,881.8 ha of additional land to the area [63]. For the Brantas/Porong delta region ( Figure 6B), we found an average shoreline change of +13.28 m/year from 2000 to 2019. The rate of shoreline change ranged from −19.98 to +111.75 m/year (EPR) (Figure 8). Unlike the Bengawan Solo delta, the Brantas/Porong delta showed no erosion at rates of more than ~20 m/year. The land-cover composition in the Brantas/Porong delta is similar to that in the Bengawan Solo delta; the acquired land is used for aquaculture ponds and as mangrove conservation areas. The massive coastal sedimentation in the delta of the Brantas/Porong river is due to a mud volcano explosion in May 2006. During this event, underground mud was lifted to the surface, carried by the river, and deposited in the delta, adding around 3,881.8 ha of additional land to the area [63].

Land Use of the Accretion and Erosion Areas in the Coastal Regions in 2000 to 2019
The mapped land-use classes in the Bengawan Solo and Brantas/Porong deltas can be found in Figure 9. Most of the delta regions were used as mangrove and aquaculture areas, especially in the Brantas/Porong delta. In contrast, in the Bengawan Solo delta, small portions of the areas were used as plantation and built-up areas. Based on the landuse classification results and the shoreline changes in 2000 and 2019, the mangrove and aquaculture areas were the main affected land-use classes in those deltas.

Land Use of the Accretion and Erosion Areas in the Coastal Regions in 2000 to 2019
The mapped land-use classes in the Bengawan Solo and Brantas/Porong deltas can be found in Figure 9. Most of the delta regions were used as mangrove and aquaculture areas, especially in the Brantas/Porong delta. In contrast, in the Bengawan Solo delta, small portions of the areas were used as plantation and built-up areas. Based on the land-use classification results and the shoreline changes in 2000 and 2019, the mangrove and aquaculture areas were the main affected land-use classes in those deltas. Based on the overlay results between the shoreline changes and the land-use information in 2000 and 2019, the type of land use that was developed in the accreted areas and lost in the eroded areas is summarized in Table 3. In the Bengawan Solo delta, around 176.59 ha of aquaculture and 40.05 ha of mangrove areas were lost because of coastal erosion; however, around 471.50 ha of aquaculture and 413.44 ha of mangrove areas were developed in the accreted areas. Similar land-use composition for the eroded and accreted areas was also found in the Brantas/Porong delta. Most of the accreted areas were used for mangrove plantations, with 251.85 ha along the shoreline of the Brantas/Porong delta, and for aquaculture, with 166.97 ha. The coastal erosion in the Brantas/Porong delta was notably lower, with only 57.4 ha, which happened mostly in the aquaculture areas.

Detecting Shorelines from Satellite Imagery Data
Various methods are available for detecting and monitoring shoreline changes from satellite imagery. These methods rely on the ability of the satellite imagery to discriminate between land and water bodies, which is commonly achieved by using (1) the histogram threshold value from water indices, such as the infrared band ratio [64] and the modified NDWI and AWEI [65,66], or (2) supervised classification [67,68]. Although histogram thresholding methods are far simpler than supervised classification, the threshold must be chosen carefully to avoid further classification noise [65]. The applicability of the histogram thresholding method to composite data of annual medians, such as the remote sensing data that we used in our study, would require further assessment. The AWEI is thought to be the most consistent index for separating land and seawater compared with the modified NDWI [65]. The assessment of which method to use should be made based on the selection of the appropriate spectral indices and threshold values, since the pixel values are taken from different dates within a year of data. In addition, the performance of each index may vary depending on the type of coast [69]. A problem with the histogram thresholding method was the rationale behind the use of supervised classification as the main method for detecting shorelines in our analysis. The use of the GEE platform in our analysis enabled the access of multi-sensor medium-resolution satellite datasets, which were useful for the purpose of classifying with greater accuracy and over a wider area of analysis.
The GMO-Maxent classification results were used in this study because of their greater accuracy when compared with RF classification, especially for the years when only Landsat data were available (e.g., 2000, 2005, and 2010). The RF method was able to provide slightly better accuracies when more data were used as input for the classification, such as the Sentinel-1 radar data for the years 2015 and 2019. Although several studies, such as the work of Farda [70], have asserted the superiority of RF as compared with GMO-Maxent, RF requires fine-tuning to optimize classification, which runs the risk of reduced accuracy, as shown in the study by Shelestov et al. [71], in which the accuracy of RF was lower than that of GMO-Maxent. Nonetheless, it should be noted that the final results of GMO-Maxent still require manual refining to obtain the most accurate shoreline results possible and to ensure the correct calculation of shoreline changes from the classification results.

Shoreline and Land-Use Changes from 2000 to 2019
The shorelines of the examined portion of East Java Province displayed a high degree of dynamism between 2000 and 2019, with an overall trend of accretion along the shoreline. Examination of the data led to the identification of two major areas of shoreline changes, namely the Brantas/Porong and Bengawan Solo deltas. In the Brantas/Porong delta, accretion was dominant, which can be attributed to sediments from the mud volcano disaster in Sidoarjo, which directed sediments to the Brantas/Porong river. During the wet periods, floods and water flows gradually flushed significant amounts of mud into the Brantas/Porong delta from the mud volcano event during the dry period [72]. Kure, Winarta, Takeda, Udo, Umeda, Mano, and Tanaka [72] also found twice the normal amount of suspended matter and particulate organic carbon in the Porong river as a result of the mud volcano event. Our study showed the gradual accretion of the Brantas/Porong delta from 2000 to 2019. However, the exact influence of the mud volcano sediments on the shoreline changes still requires further assessment, as accretion and erosion of the shoreline may also be influenced by ocean wave patterns. Nonetheless, Wijayanti [73] indicated that there was a correlation between the amount of suspended matter and the coastline changes in the Brantas/Porong delta.
The second area that showed massive changes was the Bengawan Solo delta. Although accretion dominated in this area, with an average rate of change of +10.2 m/year, some parts of the delta also suffered from massive coastal erosion, with a maximum rate of −87.44 m/year. Our results complement previous analyses, which identified erosion only in the western part of the delta from 1972 to 2016, albeit at a much slower rate (−5.54 m/year) [29]. Although our study found some major coastal erosion in the Bengawan Solo delta, this was not the case during the late 1980s. Hoekstra [74] identified the Bengawan Solo and Brantas/Porong deltas as rapidly growing deltas, which received substantial quantities of sediments, transported during the wet seasons from 1890 to 1985. The Bengawan Solo delta was formed by a single river channel, with a growth of +70 m/year [75], a figure that drastically dropped to +10.2 m/year from 2000 to 2019.
Although the observed massive coastal erosion occurred in non-urban areas, the shoreline dynamics in this area evidently require proper monitoring to avoid potential hazards in the future and to best manage the accreted land. Here, the major accretion and erosion occurred in less populated regions of the deltas, reducing the risk to society; however, coastal erosion occurred in areas of aquaculture land use, which could pose an economic problem to local communities that depend on aquaculture production [76] or which may produce a conflict in the management of the accreted lands. To minimize the socioeconomic impact of coastal erosion, coastal defense mechanisms supported by integrated coastal management should be implemented. In Indonesia, afforestation by planting mangrove in the accreted areas is the common procedure for coastal land management [77]. However, controlling accretion may be more important, especially if an offshore reef ecosystem is present, since terrestrial sediments reduce the water quality and increase the turbidity, which in turn limits the resources needed by coral reefs to grow [78]. Therefore, the detection of regional-level shoreline dynamics can be utilized to identify vulnerable coral reef ecosystems and mitigate the negative impact of ongoing sedimentation.
In addition, our results can also be used to create a measure or index for quantifying the coastal vulnerability risks by employing the direction and trend of shoreline changes to identify the elements that are at risk of shoreline changes. The coastal risk index (CRI-Med), developed by Satta et al. [79], and the coastal vulnerability index (CVI), developed by Gornitz et al. [80], use sea level rise, drought, and wave height as the hazard components and tourism, population, and land use as the elements at risk. The elements at risk, i.e., land use and population, can be identified, and the vulnerability can be calculated by looking at the trend of shoreline changes and the distance of the erosion sites to the population or the land-use areas.

Conclusions
Our study demonstrated the use of cloud computing-employing the GEE, multisensor and multi-temporal satellite imagery, supervised classification, and GIS analysis, as well as DSAS-in monitoring the shoreline dynamics from 2000 to 2019 of a large part of East Java Province, Indonesia. In our study, we found that GMO-Maxent methods can give relatively stable, high accuracies (84% to 95.2%) when compared with RF methods, especially when only Landsat data were used. The manually refined shoreline generated by the GMO-Maxent method showed general accretion, at rates of +4.12 (EPR) and +4.26 m/year (WLR). Massive changes were found in two deltas: the Brantas/Porong delta, which underwent major accretion, and the Bengawan Solo delta, which experienced both major accretion and major erosion. Most of the land areas derived from shoreline changes, such as accretion zones, were converted into mangrove and aquaculture areas, replacing the same land-use type that disappeared because of coastal erosion. The shoreline dynamics in this area require proper monitoring to ensure the maximal environmental benefit from the accreted areas and to avoid the risk of environmental hazards and socioeconomic problems due to coastal erosion.