Eleven years of mangrove–mudflat dynamics on the mud volcano-induced prograding delta in east java, indonesia Integrating uav and satellite imagery

: This article presents a novel approach to explore mangrove dynamics on a prograding delta by integrating unmanned aerial vehicle (UAV) and satellite imagery. The Porong Delta in Indonesia has a unique geographical setting with rapid delta development and expansion of the mangrove belt. This is due to an unprecedented mud load from the LUSI mud volcanic eruption. The mangrove dynamics analysis combines UAV-based Structure from Motion (SfM) photogrammetry and 11 years (2009–2019) satellite imagery cloud computing analysis by Google Earth Engine (GEE). Our analysis shows unique, high-spatiotemporal-resolution mangrove extent maps. The SfM photogrammetry analysis leads to a 3D representation of the mangrove canopy and an estimate of mangrove biophysical properties with accurate height and individual position of the mangroves stand. GEE derived vegetation indices resulted in high (three-monthly) resolution mangrove coverage dynamics over 11 years (2009–2019), yielding a value of more than 98% for the overall, producer and consumer accuracy. Combining the satellite-derived age maps and the UAV-derived spatial tree structure allowed us to monitor the mangrove dynamics on a rapidly prograding delta along with its structural attributes. This analysis is of essential value to ecologists, coastal managers, and policymakers.


Introduction
Mangroves are trees or large shrubs that grow in or adjacent to the intertidal zone and are distributed along (sub-)tropical coasts and estuaries [1,2]. Mangroves are well-known for providing a range of ecosystem services such as timber production, carbon sequestration, soil formation, nutrient cycling, habitat creation for marine and terrestrial species, and protecting coastlines by attenuating waves and limiting erosion [2][3][4]. However, land use conversion has caused mangrove forests to decline at a rate three to five times larger than the average forest loss. Hence, mangroves' important ecosystem services will be diminished [2].
In order to gain a better understanding of the value of important ecological services provided by mangroves, further efficient collection of data, e.g., mangrove extent, height, individual position, and species membership are necessary. Earth observation offers a method for the large-scale monitoring and assessment of the environment, especially when it comes to the mangrove forest inventory and spatial extent monitoring that are generally Figure 1. The study area is described in a sequence: (A) Indonesian border (light green) with East Java Province depicted in dark green, (B) East Java Province and the capital city Surabaya (represented as a red circle), (C) LUSI (lumpur (mud) and sidoarjo (the regency name)) mud volcano represented as a black polygon, and the Porong River as the black line flowing from the west to the east part of the map, and (D) Porong Estuary with LUSI island and the delta lobes.
The Brantas river originates from the volcanic complex of the Semeru and Arjuno Mountains [32]. Erosion rates on the slopes of Mount Semeru are among the highest recorded in the world (105-106 m 3 km −2 a −1 ) [33]. Sediment yield in several drainage basins of Mount Semeru is dominated by rain-triggered events during the wet season [30,33]. The Brantas watershed is densely populated and affected by anthropogenic activities such as deforestation, intensive agriculture (mainly rice cultivation), and industries [34]. Due to the geological conditions consisting of the presence of easily erodible soils and high anthropogenic activities, surface erosion is high [31]. The Porong River drains off high sediment loads, causing a prograding delta [34,35] with a progradation rates of approximately 0.4 × 10 6 m 2 y −1 over the 1935-1981 period [21,31].
In addition to the sediment load due to runoff from the hinterland, the Porong River has been experiencing an extreme sediment load due to the mud volcanic eruption in Sidoarjo, Indonesia. The LUSI mudflow is reported to be the 'largest mud eruption in the world' [36], about 18 km west of Porong Delta (Figure 1). On 29 May 2006, the boiling mud erupted at a peak flow rate of up to 180,000 m 3 day −1 [22] which declined to 50,000 m 3 d −1 in September 2011 [23]. Sixty thousand residents were forced to evacuate, and 7 km 2 of residential area was submerged with mud [36] (Figure 2a). The excessive LUSI is still actively erupting material, gas, water, clasts, and oil, albeit at a considerably reduced rate [36]. The continuous discharge of mud has been diverted to the Porong River since 2007 [23]. To reduce damage to the nearby community and environment, the mud in LUSI is first stored in a reservoir contained by 10-m-tall dyke and then diluted and disposed of by pumping to the river [21,36]. This operation has increased sediment concentration and loads of the Porong River by a factor of 3-4 compared to pre-LUSI conditions [21,23,34]. As shown in Figure 2b, due to the mudflow, the delta is rapidly prograding along with the development of mangrove belts. To date, the LUSI mud volcano is still erupting with no end in sight [36]. The delta lobes ( Figure 1) have naturally developed, while LUSI Island was created as a spoil bank to contain sediment dredged from Porong River, after severe siltation due to the mud diversion operations. It was constructed between February and November 2009 by building a 4 km series of geotube breakwaters which were attached to the natural existing Sarinah Island [23]. In 2009, 5000 Avicennia spp. seedlings were planted at this newly created wetland [23,37], and the Ministry of Marine Affairs and Fisheries continued planting thousands of Rhizophora spp. seedlings between 2010 and 2011 [26,27].

UAV Data Collection and Processing
The fieldwork was conducted at the end of the dry season (October-November) in 2019, focusing on two delta lobes (northern and southern deltas) with a total area of approximately 0.3 km 2 . This timeframe was expected to have a higher probability of clear satellite images with limited cloud cover, and more sunlight as well as less shadows to reach an optimum condition for unmanned aerial vehicle (UAV) image acquisition. UAV/ drones are widely accepted as standard survey tools in many environment settings [38], generating high-resolution data in a safe, straightforward, and cost-effective way [39]. Coastal environments are challenging for the UAV-based surveying method because of the low texture and contrast of the bed surface [39]. Careful planning of the GCPs placement and flight path can lead to reach centimetres of vertical accuracy [15,39].

Data Acquisition
Several studies demonstrated that consumer-grade UAV could achieve vertical accuracies in the order of a few centimetres to a few decimetres in coastal topographical surveying [38][39][40]. In comparison with other low-cost platforms (e.g., kite, pole, and fixedwing), a consumer (rotary) drone system with its integrated positioning system, inertial measurement unit, and stabilised camera offers flexibility and efficiency while attaining the accuracy in coastal area application [39,41,42].
In this study, we used consumer-grade UAV DJI Mavic Pro (DJI, Shenzen, China) during the field campaign. The drone has four propellers and a built-in true colour camera. The camera is equipped with 1/2.3 CMOS sensor with total effective pixels of 12.35 M, which produces a 4000 × 3000 image resolution and equipped with an electronic shutter. Information on shutter type is important since it will affect the setting on the camera calibration to compensate the rolling shutter issues as in the electronic shutter-type cameras. The overall flight time in optimal condition is 21 min with 15% remaining battery level [43].
The DroneDeploy web app (DroneDeploy, San Fransisco, CA, USA) was used to define the flight path [40]. It was planned with a flight altitude of 60 m, an overlay of 80% front overlap and 75% side overlap. Each flight was designed to cover 0.02 km 2 area with 15 min flight time to limit one battery per flight. An enhanced 3D mode was activated to improve 3D structures quality which will capture an oblique image of the objects, facing toward the inner centre of the target by carefully not to include horizon in the shots [44]. We added a buffer zone approximately 20 m apart from the edge of the low tide limit that was cropped during the processing later on to prevent interference with the SfM Photogrammetry processing [15]. With these settings, we were able to build 15 grids covering around 0.3 km 2 . All the grids created in the web app then can be synced into the mobile app to manage the flight and photo acquisition on-field automatically. The flights were conducted between 09:00 and 12:00 to get optimum natural sunlight, primarily limit the appearance of shadows on the photos. Additionally, this short time window was chosen to avoid high variation in sun intensity [45] and it covered the low-tide period. Before the flight, we placed the custom made printed red cross-shape tarpaulin as a ground control point (GCP). The GCP size is 1 m × 1 m rectangle in the inner side and 1 m diagonal cross in the outer side to make it identifiable from the 60 m altitude. We placed GCPs over the mudflat and in the middle of the delta, which has low vegetation density and registered them with DGPS in RTK mode.
In total, the properties and location of 69 mangrove trees were measured concurrently with the mudflat topography. We recorded diameter at breast height (dbh), height, and location for each tree. The dbh was recorded using a measuring tape at 130 cm above the ground. Tree height was recorded with laser rangefinder by measuring horizontal distance and hypotenuse from surveyor location to the tree stem and treetop, respectively. For the trees that were located in the proximity of the mangrove edge we recorded the position with DGPS. Trees that were located more in the centre of the forest were recorded with GPS (horizontal accuracy ± 3.5 m) because it is difficult to walk inside the forest with heavy equipment. Mudflat topography was measured with DGPS in RTK mode. The mangrove trees dataset was used as control points and groundtruthing for the created Digital Surface Model (DSM).

Data Processing
We estimated the tree's location and height based on point clouds created by the Structure from Motion (SfM) photogrammetry method. Three-dimensional point clouds derived from SfM photogrammetry produced conservative and realistic measures of tree heights [12]. Alongside with point clouds, a DSM, Digital Terrain Model (DTM), and an orthomosaic were also generated from these processes [13].
The workflow comprises of three phases, namely pre-processing, processing, and post-processing. The workflow in point clouds generation and processing is illustrated in Figure 3. Pre-processing is mainly related to data acquisition including flight route planning, GCP placement and processing of overlapped images with SfM Photogrammetry. After the data is acquired, we process the images in commercial SfM photogrammetry software Agisoft Metashape Professional 1.6.1 (Agisoft LLC, St. Petersburg, Russia). The workflow is in general as follows: (1) photo alignment, (2) photo marking, (3) dense point clouds generation, (4) exporting dense point clouds, and (5) DSM and orthomosaic generation. A detailed description of this workflow can be read in the Appendix A, Table A1.
Since the UAV camera cannot penetrate the mangrove canopy, the DSM includes vegetation and any above ground covers. Therefore, on the next workflow, we processed the point clouds to estimate the DTM. Orthomosaic was created by using surface information provided by the DSM and orthorectifying of the overlapping images.
The processing phase consists of nine steps, which are quality checking, clipping, indexing, tiling, sorting, noise removal, ground classification, height normalisation, and generating the Canopy Height Model (CHM). Details in this phase are described in Figure 4. We used LASTools (rapidlasso GmbH, Gilching, Germany) to clean the artefacts. The Cloth Simulation Filter (CSF) algorithm [46], which is efficient to extract the bare earth in lidR package [47], was employed for ground classification. As the first step, we assess the quality and information, e.g., number of points, projection, and point density with lasinfo module. Dense point clouds derived from SfM photogrammetry still contain noise and uncertainties [48,49]. Therefore, to get an optimal application in the mangroves environment, a cleaning procedure should be considered. The waterbody is the main noise contributor of the SfM photogrammetry point clouds. We visually inspected the orthomosaic and draw a polygon which excludes water and selected the area of interest. The raw point clouds were then clipped based on the polygon with lasclip module. Next, we tiled and indexed the clipped point clouds into a smaller tile of 40 × 40 m and added a buffer that was 20% of the tile size. The tiling procedure is useful to decrease the computational time by taking advantage of parallel computation. Another consideration is that, the free licensed LASTools has a point limitation about 1~15 million points [50]. Therefore, tiles should be adjusted with the allowable number of points as attached to the license. The tiled points were then sorted to rearrange the points into a space-filling curve order. Afterward, we reclassified the highest points, find the highly isolated points in the dense forest, to create a temporary ground classification. Next, we mapped all the points located 0.2 m below the temporary ground as noise. These high and low noise were removed with the lasthin and lasnoise module. The cleaned photogrammetry point clouds then were processed with the CSF algorithm in the lidR package to define the ground points. An exhaustive manual revision of the ground-no ground classification was done in the dense forest to avoid ground misclassification. The point clouds were created only from what is visible to the camera within the path not penetrating the canopies, because of that only vegetation tops/canopies were included. Therefore, careful manual revision had to be made. Finally, points were height-normalised by replacing the height of each point to the relative height of the ground-classified points. Subsequently, we converted the height-normalised points into CHM and the ground-classified points into DTM. In the post-processing phase, we used a CHM which is a gridded canopy model with 5.3 cm resolution. To locate individual trees in this very high-resolution CHM, the lidR 'tree_detection' provides two algorithms-Layer Stacking [51] and Local Maximum Filter (LMF) [52]. We set the LMF algorithm parameters window size to 5 m [15], a minimum height of 1.37 m, and a circular moving window shape which represent the crown of a mangrove tree. The resulting data CHM were tree location (x,y) and tree height (z). The detected trees were manually checked to avoid errors [12].

Tree Detection Validation
The latest forest inventory was conducted by the Ministry of Marine Affairs and Fisheries (MMAF) in 2010 [53,54]. Therefore, we could not directly validate presence of individual trees detected in CHM. Most of the tree inventory collected during the fieldwork was located around the edge of mangrove forests. Therefore, we performed visual interpretation undertaken by the three researchers [15]. We pre-selected 30 × 30 m plots reflecting the plot size as in traditional forest inventories. For each of the two delta lobes three plots were chosen, one being located at the centre, another at the northern and still another at the southern edge of the forest. The background of the researchers is (1) a groundwater engineer with no experience in ecology and remote sensing, (2) a hydraulic engineer with basic experience in remote sensing and no knowledge in ecology, and (3) a coastal engineer with basic knowledge in ecology and advanced knowledge in remote sensing. The horizontal accuracy of the detected trees and validation data can be quantified with widely used guidelines by the National Standard for Spatial Data Accuracy (NSSDA) [55] to measure the positional accuracy of the spatial datasets [56].

Satellite Data and Processing
Google Earth Engine (GEE) has been used as a tool to perform climate and hydrology and natural disaster analysis, and image processing, land use/land cover classification (LU/ LC), and urban planning [57,58]. GEE consists of petabytes of science-ready datasets and is equipped with high performance computing that can be accessed through an application programming interface (API) that is available in JavaScript and Python [59]. Some notable applications of GEE are, for instance, the Global Mangrove Watch [10], global forest cover change [60], global surface water changes [61,62], global shoreline changes [63], and continent level agricultural mapping [64]. To understand the mangrove extent development, we processed satellite imagery, such as Landsat 7/8 and Sentinel 1/2, in GEE with Python package geemap [65].

Available Dataset
We used combinations of satellite imagery from four satellites constellations available and science-ready in the GEE platform, i.e., Landsat 7 (L7), Landsat 8 (L8), Sentinel 1 (S1), and Sentinel 2 (S2). S1 and S2 are satellite constellations developed by ESA (European Space Agency) along with other constellations such as Sentinel 3 and 5. S1 and S2 satellites hold the sensors that are suitable for the mangrove classification study. S1 satellite carries a dual-polarisation C-band Synthetic Aperture Radar (SAR) instrument [66,67]. Each S1 scene in GEE has been pre-processed with the Sentinel-1 Toolbox and is science ready [68]. We employed Sentinel-1 SAR Ground Range Detected (GRD) products with 10 m spatial resolution with available dataset in GEE ('COPERNICUS/S1_GRD') from October 2014. The S2 satellite carries a Multispectral Instrument (MSI) with 13 spectral bands: 10 m resolution RGB and NIR, 20 m red edge and SWIR, and 60 m atmospheric band [69,70]. The S2 MultiSpectral Instrument (MSI) Level-1C products were used with available dataset in GEE ('COPERNICUS/S2') from June 2015. L7 and L8 are satellite constellations developed by a joint program of the USGS and the NASA. Within the GEE platform, L7 scenes have been atmospherically corrected using LEDAPS [71,72] and L8 using LasRC [73,74] (USGS Landsat Surface Reflectance Tier 1). Clouds, shadow, water, and snow were masked with the CFMASK algorithm [71,73,75]. The L8 Surface Reflectance Tier 1 ('LANDSAT/LC08/C01/T1_SR') which is available from April 2013 onwards and Landsat 7 Surface Reflectance Tier 1 ('LANDSAT/LE07/C01/T1_SR') dataset available from January 1999 onward were employed in this study.
Seasonal varying precipitation, flood control measures, and diversion operations are likely to influence the sediment yield. Mangrove flowering and seedling dispersal periods are also influenced by the season. By considering those conditions and dynamics on the delta, we proposed to map mangrove extent on a seasonal time scale. We mapped the mangrove extent that represents dry season (May and August) and wet season (November and February). A composite from one-month scenes of the particular month was created with the median value of the selected bands and indices for mangrove classification for optical sensor (S2, L7, L8). A mean function was applied for the satellite with SAR instrument (S1) [76,77]. The median value was chosen because it is less affected by outlier values that arise, for instance, from pixels affected by clouds or snow during the masking procedures [75]. To optimise the dataset availability and account for mangrove-mudflat dynamics, we employed the combination of S1 and S2 during the period of November 2015 to November 2019. The combination of optical S1 and S2 SAR was reported to improve the classification accuracy [12,76,78,79]. L8-based classification was for the period of August 2013 to August 2015, and L7-based classification was from dry season 2009 to dry season 2012. An exception has been made for L7-based classification. The L7 series of scenes has gaps due to Scan-Line Connector (SLC) failure or stripping problem [80]. Therefore, instead of using the one-month scenes we applied median values of full seasons (six months for each) during 2009-2012. Since in that period the diversion operation had just begun and the LUSI Island reclamation project had been conducted, we observed less mangrove succession in the delta. A presentation of the seasonal classification with a single map for the period of 2009-2012 is considered sufficient. All available datasets in a particular classification period were used to generate a cloud-free composite and improve the classification accuracy [76,78].

Vegetation Indices
We assessed the mangrove vegetation cover by way of four optical-related vegetation indices that are widely used in land cover characterisation [14,78,81,82], i.e., NDVI (Normalised Difference Vegetation Index) [83], NDMI (Normalised Difference Moisture Index) [15], EVI (Enhanced Vegetation Index) [84], and SAVI (Soil-Adjusted Vegetation Index) [85]. The S1 images were first pre-processed with a speckle filter (Lee refined) at a window size of 7 × 7 pixels [81,86]. A ratio channel (VV/VH) from the backscattering was generated. Lastly, the mean value of all S1 images was employed. This mean function makes the S1 composite less susceptible to variation in image acquisition [76,77]. With this additional dataset, each S1 mosaic had two bands and one index, while S2, L7, and L8 had 10, 6, and 7 spectral bands, respectively, and each was composed of four vegetation indices (Table 1). Table 2 specifies the vegetation indices formulas used in the analysis. Figure 5 describes the flowchart of the image processing procedure.

Land Cover Classification
A supervised land cover classification with machine learning algorithm was employed in GEE. Several of those algorithms have been embedded in GEE. These classifiers are Classification and Regression Tree (CART), Random Forest (RF), NaiveBayes, and Support Vector Machine (SVM) [59]. Wetland mapping, as in mangrove swamps, is one of the most challenging areas for remote sensing. RF supervised classification is reported to have the highest accuracy among widely used machine learning algorithms, e.g., K-Nearest Neighbour (KNN), SVM, Maximum Likelihood (ML), and CART [76,87]. RF comprises a collection of tree-structured classifiers to make prediction [88]. RF is more robust to noise and size reduction of the training set than CART [89], easier to implement than SVM [90], and is particularly suitable to handle high dimensional remote sensing data [91].
The land cover classification was conducted using the following workflow: (1) training data collection, (2) initiating classifier and adjusting parameters, (3) train the classifier with training dataset, (4) classify the image based on the trained dataset, and (5) accuracy checking. In the first step, the training dataset is provided from the very high-resolution orthomosaic created on the UAV data processing. We created three classifications (water, mangroves, and mudflat area) by plotting region of interests (ROIs) which represent these classifications. As the intertidal mudflat varies over time by tidal variation, proper consideration should be made by not only referencing the mudflat's ROI based on the orthomosaic but also considering the median composite in November 2019. The ROIs polygons were made in QGIS as a shapefile and uploaded to GEE as assets. Next, sample points were created based on supervised ROIs. Accuracy of RF classification is sensitive to the sample size and spatial distribution [92]. Within each stratum, which is based on supervised ROIs, the stratums are randomly sampled [93]. Stratified random sampling function in GEE was used to sample the training data. Cochran's formula (1) was employed to determine sample size by assuming an unknown proportion for each class [94].
Here, n 0 is the sample size per class, p is the proportion of the population which has the class in question, q = 1 − p, Z is the z-value for the given confidence, e is the margin of error. The samples were subdivided into a 70% training set and 30% of validation [88]. In the next step, to optimise the computational performance, RF classifier with 200 decision trees [78] was initiated and trained as several studies suggested 100-500 as the optimal number [89,91,[95][96][97]. The trained dataset then is used to classify the composite.

Accuracy Assessment and Validation
We used accuracy assessment functions embedded in GEE for those parameters. The final ground-truthing is based on the VHR orthomosaic derived from UAV image acquisitions in November 2019. We used the confusion error matrix with these parameters: Overall Accuracy (OA), Kappa Coefficient, Producer Accuracy (PA), and Consumer Accuracy (CA) [87,90]. The OA is calculated by summing the number of correctly classified values and divided by the total number of values. PA is determined by comparing the number of correctly classified values of a particular class and number of reference pixels of the same class, while CA can be calculated by dividing the number of correctly classified values of a particular class and number of classified pixels in the class [98]. The OA, PA, and CA are expressed as percentage, with 100% accuracy representing a perfect classification. Kappa measures the difference between the actual agreement in the error matrix and the chance agreement that is indicated by the row and column total [93]. A kappa value of 0 represents no agreement, and a value of 1 indicates perfect agreement.

Results
This section outlines the result of drone-based point clouds, tree detection and validation, and satellite-based mangrove extent change over time.

Point Clouds
In total, we processed 2860 UAV images (2020 images for north delta lobe and 840 images for the southern delta lobe), which is equal to an average of 200 images per grid. SfM Photogrammetry provided point clouds, DSM, and orthomosaic for each delta lobe. The generated raw point clouds had 400.5 million and 136.6 million points, which correspond to average point densities of 951.25 m −2 and 736 m −2 for the northern and southern delta lobes, respectively. The derived DSM has a resolution of 5.33 cm/pixel and an orthomosaic with 2.66 cm/pixel resolution. The products in total covered an area of 0.44 km 2 and were corrected based on GCPs with a total error of 0.06 m. The discrepancy between the planned flight area and the product is generated from other areas captured during the UAV image acquisition.
A small subset of the point clouds is shown in Figure 6. The set of figures illustrates the step-by-step point clouds processing. Figure 6a shows subtle noise yielded on the SfM-based point clouds. The raw point clouds were refined by classifying the high and low noise, resulting in cleaned points as depicted in Figure 6d. The ground-classified points with the CSF method were then evaluated, especially those situated below the dense vegetation since SfM-based point clouds provide no information below the canopy cover. Hence, an exhaustive manual revision was conducted by checking and correcting the classified ground points. The final point classification categories chosen for the analysis were non-ground (1), ground (2), high vegetation (5), and noise (7). The final result of the point cloud processing is a height-normalised point clouds that can be described with all elevations were normalised with respect to the ground, i.e., an elevation of 0 m. Figure 7 shows DTM and DSM derived from the final classification of the point clouds.

Canopy Height Model (CHM) and Tree Detection
The height normalised point clouds were rasterised to a CHM with a resolution of 0.1 m, and subjected to the individual tree detection. The 'tree_detection' function in the lidR package with the LMF algorithm was adjusted to detect mangrove trees above the breast height. The individual trees derived from the algorithm can be observed in Figure 8. The root mean square error (RMSE) analysis of the individual tree detection based on fieldwork data and visual inspection is described in detail in Appendix B, Table A2. The RMSE value for tree location was 0.23 m on average. All three validators provided similar RMSE r with values ranging 0.15, 0.28, and 0.25 m, while the maximum RMSE r was 0.48 m. As observed, the number of the detected trees for the dense mangrove forest was underestimated, whereas it was more accurate in the sparse mangrove forest.
The CHM-derived height demonstrated that the mangroves' median height is 3.5 m in the southern delta while the northern delta with its older mangroves had a median height of 4.2 m. As shown in Figure 9, the northern delta had bi-modal tree height distribution. The distribution is likely correlated to the mangrove planting activities in 2016. Corresponding rectangular shaped mangrove areas can be clearly seen in Figure 7a.     Figure 11). On the northern delta lobe, mangroves appeared in 2016 followed by the southern delta lobe in 2018. The land conversion to fish ponds in the hinterland is apparent as well in the figures. Figure 10 provides an illustrative visualisation of the mangrove dynamics. Starting in 2006 the mangrove forests tended to expand seawards. The mangroves in northern and southern delta lobes contributed significantly to this. Considering the time series map, it is likely that the southern delta mangrove will be attached to the LUSI Island within the next few years.
The time series of the mangrove extent has been extracted for both on the region of interest (ROI) (Figure 12a) as well as for the LUSI Island-delta lobes (Figure 12b). Both areas exhibit a similar positive trend of development. Generally, it varies with season. The mangroves area recedes during the transition from dry to wet season and the mangroves regrow during the wet to dry season. However, the delta lobes' mangroves have a slightly different seasonal pattern. As in Figure 12b, the mangrove's seasonal variation is less pronounced. It is likely because it is more isolated and is affected less by human activities.

Accuracy Assessment of Porong's Mangrove Classification
Four confusion matrices of L7, L8, and S1/2 based classification show high accuracy when compared with the UAV ground reference data. The reference data was subdivided into 70% fraction of the total sample points (cf. Section 2.3.3) for training purposes, equal to 808 points for all classes. The rest (30% of the sample points) or 347 points in total were used for validation samples. The values of Overall Accuracy (OA), Kappa Coefficient, Producer Accuracy (PA), and Consumer Accuracy (CA) of all the classifications were all above 98%. L7 or L8 alone are already excellent sources for the classification of the onemonth composite, but the combination of S1 and S2 is even superior to those. The trained OA, Kappa, PA, and CA of S1-S2 classification were all 100%. Due to the SLC error, the trained L7 classification generated from a six-month mean composite also resulted in 100% OA, Kappa, PA, and CA values. The L8-trained OA, Kappa, PA, and CA were 99.88%, 0.99, 99.65%, and 99.62%, respectively. The validation accuracies also showed consistently excellent values similar to those of the training. The combination of S1-S2 has a value of 100% for all confusion matrices. The L8-validated OA, Kappa, PA, and CA were 99.39%, 0.99, 97.92%, and 98.39%, respectively. The L7-validated OA, Kappa, PA, and CA were 99.54%, 0.99, 97.83%, and 98.18%, respectively.

Age Map
The age map (Figure 13) was estimated and referenced to November 2019 and derived backward to 2009. The age map indicates that the mangrove expansion of forests being attached to the mainland began in 2014. In comparison, the isolated forests (on the delta lobes) were found to have started expanding in 2016, most likely initiated by the mangrove planting in the northern delta, in contrast to the natural mangrove succession that took place in the southern delta. By taking advantage of the age map and the high-resolution CHM, a relationship of mangrove height dependent on stand age was setup for the two Porong delta lobes. A mean height distribution across the mangrove age map was calculated by clipping the mangroves' age polygon to a CHM raster. Figure 14 compares the age-height relationships from the delta lobes' mangroves. The northern delta mangroves were consistently taller than the southern ones at the same stand age. By only taking the mangrove trees into account that were taller than 1.3 m, the average annual mangrove growth of the northern delta amounted to 2.26 m yr −1 . In contrast, the mangroves in the southern delta had an annual mangrove growth rate of 1.71 m yr −1 .

Discussion
The study aimed to quantify the mangrove dynamics arising from the massive LUSI mudflow diversion operation that followed the extreme mud volcano eruption, in Sidoarjo, Indonesia. Our investigation started in January 2009 and continued until November 2019. This investigation of mangrove dynamics combined usage of UAV-based and satellite analyses with GEE cloud computing. This approach included the successful retrieval of mangrove biophysical properties in terms of canopy height and the individual position of mangrove trees (Figure 8) as well as a time series of mangrove belt development ( Figure 10). This approach resulted to a new set of mangrove extent maps for our study location that is not covered by other existing products such as the WMA [1] and GMFD v1 [8]  . Therefore, from those datasets the mangrove dynamics in Porong Estuary cannot be deduced. The spatial and temporal resolution of those dataset is not sufficient as Porong has been experiencing extreme mud influx from LUSI and its transformation has been so rapid. This study resulted in a higher temporal resolution and detailed three monthly mangrove classification and their biophysical properties, especially in the highly dynamic delta lobes that expand seaward.

UAV-Based Mangrove Forest Inventory
The derived tree locations were in close correspondence with field inventory data and expert analyses at an RMSE r in the order of 20 cm. However, since the low-cost UAV is equipped with an RGB camera system, the estimated tree location is limited to the description of the canopy top inherent in the CHM. It is different from the LiDAR system where the beam is able to penetrate the dense canopy cover. The UAV-based SfM photogrammetry method tended to underestimate the number of individual trees, even though it has a higher accuracy in the sparse mangrove forest. Based on the result, the low accuracy in the dense forest of the northern delta lobe is likely due to the mangrove plantation programme that has created a more homogenous canopy and an almost flat CHM. This prevented complete detection of all treetop positions. However, despite the underestimated tree location, considering that colonising mangrove species during primary succession is shade-intolerant, therefore, the loss in detection might be marginal, in particular when the canopy is still sparse. The estimated individual positions of the trees can give us more information regarding the height distribution of the mangrove's trees and density (Supplementary Materials S1). As shown in Figure 9 we can observe the influence of mangrove planting on the north delta to the height variation indicated by the bimodal distribution (Figure 9a).
Despite its limitation, the UAV-based tree detection offers the advantages of lower technical requirements and lower costs when compared to a LiDAR system [39]. With miniaturisation, increasing proliferation, and advancement in sensor technology, the possibility arises of adding multispectral sensors to accuracy of the UAV-based method while keeping it lightweight and low-cost. The image acquisition was conducted for five days or equal to 8.8 hectares/day since it was limited by the time frame of the optimum sunlight and battery capacity. The UAV survey in this study was conducted by two people. In comparison, traditional mangrove forest inventory with three experienced surveyors needs two complete days (seven hours/ day) to sample 0.05 hectare of the forest [15]. Thus, UAV-based forest inventory is likely to increase the possibility for frequent and rapid mangrove monitoring. Moreover, this method can be further developed as citizen science monitoring. With an off-the-shelf drone, pre-planned flight, open tutorial, and validation procedure, the public can be involved in this kind of mangrove monitoring programme. Depending on the needs, the UAV-based SfM photogrammetry is useful to estimate mangrove biophysical properties from the plot level up to the landscape-level. Moreover, another product generated in the workflow, for instance, the DTM is useful in studying the biogeomorphology to understand the feedback loop of vegetation and the environment forcing.

Mangrove Belt Expansion Identification in Google Earth Engine
The expansion of the mangrove belt has been substantiated with the supervised land cover classification. We have evidenced a total of 11 years of mangrove development (2009-2019) with high accuracy (more than 98%). Cloud computation with GEE to our advantage, we generated three monthly mangrove maps, with the exception of sixth-monthly maps during the period of 2009-2012. Most importantly, our map series has higher frequency than the above-mentioned global mangrove products.

Seasonal Pattern of Mangrove Dynamics
Based on the resulting maps, it is apparent that the year 2016 marked the start of the positive mangrove expansion, which seen in Figure 10 and the time-lapse animation in the Supplementary Animation S1. There were two of the largest expansion areas in the Porong River mouth contributing 24% of the total mangrove area extracted inside the ROI as of November 2019 ( Figure 12). It is highly likely that the mudflow diversion procedure has promoted mangrove belt expansion in Porong Estuary. Mud pumping disposal of LUSI discharged extra sediment supply thereby providing a suitable environment that enhanced mangrove expansion. Figure 11 provides detail of the significant increase of mangrove extent both on the ROI and focusing on the estuary. Dredging operation and the completion of LUSI Island construction in 2011 have contributed to the prominent increase of the mangrove area.
The mangrove expansion followed a seasonal pattern as Figure 12 clearly indicates. We observed recession of the mangrove extent during the transition of dry to wet season and regrowth during the wet to dry season. The wet season has the highest sedimentation rates related to the mudflow pumping operation in LUSI. The mud is disposed of in the wet season and stored in the reservoir during the dry season. The seasonal pattern in Porong was also indicated in the study by Sidik et al. (2016) [21]. In principle, net development trend of the mangrove area is positive. Investigated further, we can see amplitude of the high-low signal differs in the period of 2013-2017 from that in 2018-2019. In 2013-2017, we can see low mangrove extent in the wet season followed by the high regrowth in the dry season. In contrast, from 2018 to 2019, we observe only a slight decrease in the wet season and a high roughly two-fold increase in the dry season. Focusing on the Porong river mouth, as shown in Figure 12b, the seasonal fluctuation is also recognised but is less pronounced as in Figure 12a. The graphs reveal continuous expansion until August 2013, irrespective of the region, thus no recession visible in both regions. From August to November 2013, there is a sudden recession, i.e., sharp decrease of area in region A and slightly less in region B. The recession took place in the middle of the dry season and the beginning of the wet season, thus at the end of the dry season. There is recovery occurring from November 2013 until August 2014, where recovery is more pronounced for region A (February and May clear satellite images are lacking in 2014). The recovery occurred during the wet to dry season. It seems that the recession occurred at the end of the dry season while recovering and return to expansion happens in the middle of the wet season and the first half of the dry season. This 2013 recession-recovery behaviour repeats in 2014-2017, where the fluctuation is noticeable in region A and less in region B. It is likely in the beginning mangroves start growing on the newly deposited, homogeneously distributed mud since its sediment attributes are rather suitable for their growth. Mangroves have access to water during wet season as well as during dry season. However, because of their existence, further sediment is deposited at the margin of the forest. This possibly transforms them into the basin mangrove type in certain areas, more so in the landward direction, less so on the delta lobes. When the trees start generative production, i.e., propagule production, at an age of three years seedling that are known to be more sensitive to salt and drought might die under the extreme conditions at the end of the dry season in the basins. This condition has developed from the interplay of vegetation and geomorphological processes-namely, the massive sediment load. Studies reported buried pneumatophore leads to mangrove mortality, for example in Mekong Delta [100]; high sedimentation concentration reduces the oxygen level in the mangrove's root [21,101]; and an enhanced of growth of Micronesia trees that is possibly due to associated decreases in root zone salinity [102].
The delta lobes' mangroves response to high sedimentation rates is generally appear as lack of growth instead of a decrease. It is likely because the delta lobes are relatively isolated and therefore less affected by human activities, such as land conversion to fishponds that is obvious in the maps. Given the sharp decrease-increase pattern in Figure 12a (period 2014-2018) and slight decrease-increase as Figure 12a (period 2018-2019) Figure 12b it is likely that anthropogenic activities played an additional role. However, as a verification of that role is beyond the scope of this article, we suggest further investigation of this topic.

Mangroves' Age Class Estimation
The age map ( Figure 13) was derived to understand the mangroves succession. When it is combined with the CHM, we can estimate the annual growth rate of the mangrove trees ( Figure 14). The age map indicates a difference in the trend of mangrove expansion for those located in the river mouth and on the mainland. Expansion in the river mouth was altered by the LUSI Island's construction and the built-up of the delta lobes. It is likely that the mangrove planting programme affected the northern delta's spatial tree height variation. Apparently, it is also visible in the age-height relationship displayed in Figure 14. As the mangroves on the northern delta were planted, the mangroves tended to be in average two times taller in their first year than its southern counterparts. However, the data used to derive the age-height graph in this article is probably not long enough; for instance, the southern delta provides only three years of observation. Thus, two growth rates for the relationship, which is insufficient for performing statistical tests. Nevertheless, mangrove growth rates had a similar trend for both locations, despite a difference likely due to the mangrove planting. Moreover, the graph demonstrates the advantage of combining UAVbased VHR data and satellite imagery to characterise the mangrove structural attributes that need a large workforce if done traditionally.

Implications of the Study
To our knowledge, this study presents the first attempt to explore the mangrove dynamics in a prograding delta setting by integrating UAV and multiple sources of satellite imagery (L7, L8, S1, and S2) in GEE. Studies to classify mangrove forest have been conducted in the last two decades [8,10]. They were based not only on one source but also on the combination of multiple satellite sources [12,14,103,104] and became more popular with the advent of GEE. The resulting three-monthly classification maps are more frequent than the commonly produced annual mangrove maps. The frequent monitoring is deemed necessary in this environmental setting, since the sediment influx continues at an unprecedented rate, promoting rapid delta development and mangrove belt expansion. The results of the UAV and the entire methodology presented here shows the advantage of using an off-the-shelf drone. The CHM and individual tree locations present important structural attributes to characterise the mangrove forest. It is essential as there is likely a much closer linkage between diameter at breast height (dbh) and height rather than between crown radii and dbh. Accurate information of mangrove biophysical structure is critical for the ecologist, coastal manager, or policymaker.

Conclusions
The LUSI mud volcano eruption, arguably the largest mud eruption in the world, certainly affects the downstream landscape development. Particularly with the diversion operation that conveys a large amount of mudflow sediment into the river that eventually stimulates a rapid progradation of the delta compared to pre-LUSI conditions. This offers a unique opportunity to analyse mangrove dynamics on a rapidly prograding delta. We created time series of mangrove extent maps and mangrove biophysical structure with the following steps: (a) biophysical properties retrieval using UAV-based SfM photogrammetry; (b) three-monthly classification of mangrove areas using L7, L8, S1, and S2 in Google Earth Engine; and (c) derivation of mangrove age maps based on the satellite imagery. This improvement enables us to capture the highly dynamic setting in the study area. Moreover, the off-the-shelf UAV offers an efficient yet accurate technique to retrieve the important structural attributes, such as individual tree location and canopy height. When combined with satellite imagery analyses, the information can be used to characterise the mangrove forest and assess the effect of excessive mudflow discharge on the delta and mangrove development.
The proposed approaches allowed us to monitor the dynamics of mangrove extent and structural attributes in a rapidly prograding delta. The random forest supervised classification demonstrated a high accuracy (OA > 99.39%, kappa value 0.99, PA > 97.83%, and CA > 98.18%). The highest accuracy was obtained in the classification from the combined Sentinel 1 and 2, while the lowest resulted from the Landsat 7. The 11 years of mangrove extent mapping provided evidence of an overall positive trend in mangrove extent overlaid by seasonal variation. The receding mangrove area is detected during the transition from dry to wet season and regrow during the wet to dry season. The individual trees height and position derived from UAV showed the different distribution of the north and south delta lobes that is likely related to the mangrove planting.
The method enabled us to retrieve rapid and accurate information on mangrove biophysical properties with an off-the-shelf drone. Moreover, in combination with GEE based cloud computing, it is possible to derive a high spatiotemporal resolution mangrove extent map. The combination of UAV-derived spatial tree structure and the satellite-derived maps are needed to support fast and frequent mangrove monitoring that will be valuable for the ecologists, coastal managers, or policymakers.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
Here we explain the settings of the pre-processing phase to generate 3D point clouds in the Agisoft Metashape Profesional. The settings in this appendix are for Agisoft Metashape Professional version 1.6.3 build 10723 (64 bit). The detailed description is depicted in Table A1 below.

Appendix B
The calculation of positional accuracy of the detected trees is following the guidelines by the NSSDA Part 3: National Standard for Spatial Data Accuracy [55]. We assessed the horizontal accuracy of the detected trees against the visual observation by the three observer and fieldwork measurement. The RMSE values were assessed using the following equations: