*Article* **Mangrove Ecosystem Mapping Using Sentinel-1 and Sentinel-2 Satellite Images and Random Forest Algorithm in Google Earth Engine**

**Arsalan Ghorbanian 1, Soheil Zaghian 1, Reza Mohammadi Asiyabi 2, Meisam Amani 3, Ali Mohammadzadeh <sup>1</sup> and Sadegh Jamali 4,\***


**Abstract:** Mangroves are among the most productive ecosystems in existence, with many ecological benefits. Therefore, generating accurate thematic maps from mangrove ecosystems is crucial for protecting, conserving, and reforestation planning for these valuable natural resources. In this paper, Sentinel-1 and Sentinel-2 satellite images were used in synergy to produce a detailed mangrove ecosystem map of the Hara protected area, Qeshm, Iran, at 10 m spatial resolution within the Google Earth Engine (GEE) cloud computing platform. In this regard, 86 Sentinel-1 and 41 Sentinel-2 data, acquired in 2019, were employed to generate seasonal optical and synthetic aperture radar (SAR) features. Afterward, seasonal features were inserted into a pixel-based random forest (RF) classifier, resulting in an accurate mangrove ecosystem map with average overall accuracy (OA) and Kappa coefficient (KC) of 93.23% and 0.92, respectively, wherein all classes (except aerial roots) achieved high producer and user accuracies of over 90%. Furthermore, comprehensive quantitative and qualitative assessments were performed to investigate the robustness of the proposed approach, and the accurate and stable results achieved through cross-validation and consistency checks confirmed its robustness and applicability. It was revealed that seasonal features and the integration of multi-source remote sensing data contributed towards obtaining a more reliable mangrove ecosystem map. The proposed approach relies on a straightforward yet effective workflow for mangrove ecosystem mapping, with a high rate of automation that can be easily implemented for frequent and precise mapping in other parts of the world. Overall, the proposed workflow can further improve the conservation and sustainable management of these valuable natural resources.

**Keywords:** mangrove ecosystem; random forest (RF); Google Earth Engine (GEE); Sentinel; synthetic aperture radar (SAR); optical; aerial roots

#### **1. Introduction**

Mangroves are unique ecosystems that grow along tropical and sub-tropical coastlines. They provide many ecological benefits, including coastal protection, carbon sequestration, and waste and pollution assimilation [1–5]. Despite their significant environmental services, mangroves continue to disappear due to anthropogenic activities and climate change [6,7]. For instance, over the last five decades, approximately 20–30% of global mangroves have disappeared due to various phenomena, such as urban expansion, conversion to aquaculture, sea-level rise, and sediment alterations [8–13]. Therefore, accurate spatial

**Citation:** Ghorbanian, A.; Zaghian, S.; Asiyabi, R.M.; Amani, M.; Mohammadzadeh, A.; Jamali, S. Mangrove Ecosystem Mapping Using Sentinel-1 and Sentinel-2 Satellite Images and Random Forest Algorithm in Google Earth Engine. *Remote Sens.* **2021**, *13*, 2565. https:// doi.org/10.3390/rs13132565

Academic Editor: Chandra Giri

Received: 19 May 2021 Accepted: 25 June 2021 Published: 30 June 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

and temporal mapping and monitoring of mangrove ecosystems are crucial for natural resource conservation and sustainable development goals [14].

In situ observations provide the most accurate information about mangroves. However, collecting in situ observations through field surveys is challenging, due to the limited accessibility of mangrove communities, as they are located in harsh and tidally inundated environments [15]. Therefore, remote sensing has been recognized as an efficient and cost-effective approach for mapping mangroves [16].

In this regard, multi-spectral [17–19], hyperspectral [20], light detection and ranging (LIDAR) [21], and synthetic aperture radar (SAR) [22] datasets have been utilized for mangrove studies. For instance, Manna and Raychaudhuri [23] examined the potential of Sentinel-2 images for mangrove mapping in Sunderban, India. Moreover, Zhu et al. [21] employed unmanned aerial vehicle (UAV) optical and LIDAR data to map mangroveinundation patterns in Fujian, China. Furthermore, Kabiri [24] utilized RGB images, acquired by UAV, to classify the coastal ecosystem of Nayband Bay, Iran, into five classes of mangrove, shallow water, deep water, vegetation, and sand. The ortho-images and reference samples were fed to a maximum likelihood classifier to fulfill this task, producing a land cover map with an OA of 87.6%. Likewise, Toosi et al. [25] investigated the suitability of combining Sentinel-2 and WordView-2 images to produce a mangrove ecosystem map with eight classes. An upscaling approach in three stages was applied to produce a wall-towall land cover map based on a single-date Sentinel-2 image, resulting in the OA and Kappa coefficient (KC) of 65.5% and 0.63, respectively. The incorporation of single-date satellite imagery along with few reference samples were two limiting factors of their research, leading to moderate accuracy.

Most conducted studies on mangroves using remote sensing data have implemented traditional approaches on local computers, requiring manual acquiring, correcting, and processing of satellite images. Therefore, Cárdenas et al. [26] encouraged scholars to employ cloud-computing platforms, such as Google Earth Engine (GEE), making mangrove mapping more efficient [27,28]. This platform enables the automation of repetitive tasks (e.g., image acquisition, calibrating, and processing) and, thus, decreasing the dedicated time up to 60% [26].

GEE is a cloud platform enabling high-performance computing capabilities for geospatial data processing [29,30]. GEE hosts petabytes of satellite images, which have been efficiently employed in a variety of Earth science-related studies, such as land cover mapping [31–33], monitoring of volcanic thermal anomalies [34], monitoring of forest health [35,36], and wildfire damage assessment [37]. This cloud platform was also used in several mangrove studies, especially mangrove extent mapping [8,15,38–40]. For example, Mondal et al. [39] combined annual downscaled Sentinel-2 images and two machine learning algorithms of random forest (RF) and classification and regression trees (CART) within GEE for mangrove extent mapping along the coasts of Senegal and Gambia. The final mangrove extent maps of RF and CART had average OAs of about 93.44% and 92.18%, respectively. Furthermore, Sentinel-1 and Sentinel-2 data were employed in conjunction to produce a 10 m-resolution mangrove extent map of China [15]. In this regard, quantile synthesis features derived from both datasets were applied to an RF classifier using both pixel-based and object-based approaches. Finally, the classification results were improved by the constraint of tidal flats and visual manipulations. It was reported that the pixel-based approach obtained a higher OA, of approximately 95%, based on two-class (i.e., mangroves and non-mangroves) mapping.

Almost all of the mangrove-related studies within GEE, to the best of our knowledge, were dedicated to mangrove extent mapping [8,15,39–41] and, thus, the potential of this cloud platform for detailed mangrove ecosystem mapping was not fully explored. Additionally, other studies that were carried out to map a detailed mangrove ecosystem implemented traditional approaches on local computers using a few satellite observations [23,25]. Therefore, in this study, the GEE cloud computing platform was integrated with open-access satellite images to produce a detailed mangrove ecosystem map, advancing conservation and sustainable management of these valuable ecosystems through a higher automation level. The proposed approach uses an uncomplicated yet effective framework, allowing a consistent mangrove ecosystem mapping workflow for monitoring and assessment purposes. In this regard, the objectives of the present study are summarized as (1) incorporating multi-source images of Sentinel-1 and Sentinel-2 for accurate mangrove ecosystem mapping within GEE through an efficient and reproducible workflow, and (2) employing dense seasonal time-series observations to alleviate the tidal effect for more accurate mapping without any further tidal refinements. The robustness of the proposed method was then examined through cross-validation and consistency analyses. Furthermore, the contribution of multi-source remote sensing data and seasonal observations were also investigated.

#### **2. Materials and Methods**

#### *2.1. Study Area*

The study area covers the mangrove ecosystem in the Hara protected area of Qeshm Island, southern Iran (see Figure 1), which is under protection by different international conventions [42]. It is approximately centered at latitude and longitude of 26◦ 50 N and 55◦ 44 E, respectively, between the northwest estuaries of Qeshm Island and Hormozgan province. This area is the largest mangrove ecosystem in the Persian Gulf and Oman Sea coasts [43]. The study area includes grooved tidal channels, wherein tides are semi-diurnal. Furthermore, this region is affected by considerable tidal fluctuations, which necessitate taking the tidal effect into account for accurate mangrove mapping. There are generally two mangrove species, *Avicennia marina* and *Rhizophora mucronata*, in the Persian Gulf, the dominant of which is *Avicennia marina*, which grows in oxygen-poor sediments [44,45]. The local community uses this mangrove ecosystem for fishing, leaf-cutting, and regular boat journeys (i.e., tourism), which have negatively impacted the ecosystem, suggesting the necessity of frequent monitoring for conservation and natural resource management. Additionally, due to its proximity to the Strait of Hormuz, through which a large number of oil tankers pass, this area is also adversely impacted by oil leakage [46].

**Figure 1.** The geographical extent of the Hara protected area of Qeshm Island in southern Iran, along with the spatial distribution of reference samples over the study area (reference samples information is provided in Table 1).


**Table 1.** The number and area of training and test polygons for nine classes that were considered in this study.

#### *2.2. Datasets*

In this section, the datasets employed for mangrove ecosystem mapping are explained. First, a description of the collection and preparation of reference samples is provided, and then the Sentinel-1 and Sentinel-2 satellite datasets are described.

#### 2.2.1. Reference Samples

In this study, precise visual interpretation depended on collecting reference samples from high-resolution satellite images available in ArcMap and Google Earth. Additionally, false-color composite satellite imagery and previous mangrove ecosystem maps were used. Homogenous sites were considered for reference sample collection to mitigate the challenge of mixed pixels by avoiding fragmented areas. In total, nine classes with adequate (i.e., in terms of land cover portion in the study area and possible complexity) reference samples and appropriate spatial distribution (i.e., distributed over the study area) were generated (see Figure 1). Reference samples were then randomly split into two groups of training (50%) and test (50%) samples. Random splitting leads to low bias in the performance of the final classification results [47]. However, the primary challenge of random sampling is the information leak between training and test samples [48]. In other words, random sampling at pixel unit causes the training and test datasets to include reference samples from the same polygons. This issue increases the spatial autocorrelation between training and test datasets, which affects the accuracy assessment results and decreases the generality of the classifier [49]. Therefore, to avoid this, the random splitting step was conducted at the polygon unit, which also spatially disjointed the training and test samples. It should be mentioned that the random splitting step was implemented ten times to enable applying a cross-validation procedure for performance evaluation, which can also prove the applicability and robustness of the proposed method for accurate and detailed mangrove ecosystem mapping [39]. Table 1 provides the number of training and test polygons (i.e., the average value in ten iterations) and their corresponding area. In total, 265 and 264 training and test polygons with an area of about 125.22 ha and 130.38 ha were generated, respectively.

#### 2.2.2. Satellite Images

The time-series Sentinel-1 and Sentinel-2 satellite images were integrated to produce an accurate mangrove ecosystem map. Combining SAR and optical data allows the detection of different physical and spectral characteristics of land covers and, thus, their integration may achieve precise classification results [50–52]. Additionally, time-series data enables consideration of the water level fluctuations and tidal effects in the mangrove ecosystem, which can also increase the reliability of the classification results [53].

Sentinel-1 is a European SAR satellite, which acquires C-band data in dual-polarization in all-weather conditions with a 6-day temporal resolution. Level-1C ground range detected (GRD) images with 10 m spatial resolution in both ascending and descending modes were employed [54,55]. In total, 86 Sentinel-1 scenes in the VV (vertical transmittance and receiving) and VH (vertical transmittance and horizontal receiving) polarizations, acquired in 2019 (i.e., from 1 January 2019 to 1 January 2020), were employed (see Table 2).

Sentinel-2 is also a European platform launched by the European Space Agency (ESA) and carries the MultiSpectral Instrument (MSI) sensor [54]. This sensor records the Earth's surface information in 13 spectral bands from visible to shortwave infrared (SWIR) regions, with different spatial resolutions ranging between 10 m and 60 m. In this study, only four bands of blue, green, red, and near infrared (NIR), which are captured with 10 m spatial resolution, were used. As it is acknowledged that satellite imagery with higher spatial resolutions improves mangrove ecosystem mapping, [25], here, only bands with 10 m spatial resolution were employed. This is mainly rooted in the fact that higher spatial resolution imagery improves the delineation of mangrove ecosystem classes, especially identifying mangrove patches with a small area or narrow shapes [56]. In total, 41 Sentinel-2 images, acquired in 2019 (i.e., from 1 January 2019 to 1 January 2020), were considered for mangrove ecosystem mapping (see Table 2).


**Table 2.** The number of Sentinel-1 and Sentinel-2 satellite images in each season.

#### **3. Methodology**

The schematic framework, providing an overview of the proposed approach, is presented in Figure 2. This section comprises three sub-sections, in which the satellite data preprocessing, classification workflow, and accuracy assessment procedure are explained in detail.

**Figure 2.** Framework of the proposed method for mangrove ecosystem mapping using the random forest (RF) algorithm and a combination of the Sentinel-1 and Sentinel-2 satellite images within the Google Earth Engine (GEE) platform.

#### *3.1. Satellite Data Preprocessing*

Sentinel-1 GRD data are available within GEE with the snippet of (Image Collection ID: *COPERNUCUS/S1\_GRD*). They are generally ready-to-use data because several preprocessing steps are initially applied to them by the GEE developers. These data are already converted to the backscattering coefficient (σ◦, dB) and are ortho-rectified. In particular, five preprocessing steps (1) applying orbit file correction, (2) GRD border noise removal, (3) thermal noise removal, (4) radiometric calibration, and (5) terrain correction were applied to each Sentinel-1 scene, while their detailed information is provided by the GEE developers (https://developers.google.com/earth-engine/guides/sentinel1 accessed on 29 June 2021). Afterward, all available Sentinel-1 scenes in 2019 were categorized based on the seasons of acquisition, followed by applying a *mean reducer* to generate seasonal time-series data. Downscaling time-series Sentinel-1 data by the *mean reducer* function produced seasonal datasets, which are less susceptible to image acquisition conditions, and reduced speckle noise [31].

Sentinel-2 top of atmosphere (TOA) reflectance data, which are available within GEE by the snippet (Image Collection ID: *COPERNICUS/S2*), were also used in this study. The TOA reflectance values were derived through radiometric calibration of raw data. Because of the importance of applying cloud masking, a filtering step was first implemented to remove Sentinel-2 scenes with a cloud cover percentage of higher than 5%. Subsequently, similar to Sentinel-1 data, a seasonal *median reducer* was applied to all remaining Sentinel-2 scenes, generating seasonal optical features for classification tasks. The *median reducer* function allows the production of cloud-free seasonal datasets, in which the noisy, very dark, and very bright pixels are also removed [31,57].

Ultimately, eight SAR features (4 VV + 4 VH) and 16 optical features (4 blue + 4 green + 4 red + 4 NIR bands) were used in synergy to produce mangrove ecosystem maps. It is well acknowledged that the quality of the classification results directly depends on the input features [58,59]. As such, the integration of multi-source (i.e., SAR + optical) data can increase the discriminative capability of the classifier [53]. Moreover, time-series satellite data can manifest the water level fluctuations in estuaries, such as mangrove ecosystems [53]. Consequently, seasonal datasets can mitigate the tidal effects in the study area and allow producing cloud-free mosaics.

#### *3.2. Classification*

Different classification algorithms have been employed for mangrove mapping using satellite images [23,39,60]. In this regard, choosing the most appropriate classifier, in addition to selecting discriminative features, is important, and directly affects the classification results. Among classifiers, random forest (RF) proved to be an efficient algorithm in mangrove mapping studies [15,44,61]. For instance, Toosi et al. [44] compared four frequently used non-parametric classifiers (i.e., RF, support vector machine (SVM) with linear and radial basis function kernels, and regularized discriminant analysis) for mangrove ecosystem mapping, and concluded that the RF classifier was superior.

RF is an authoritative non-parametric classifier, which employs the bootstrap aggregation technique to combine the classification results of various independent random decision trees and to predict the class label [62]. Each of these random decision trees is trained by a subset of training samples, called in-bag samples, and uses the remnant, called out-of-bag samples, for internal cross-validation. Later, their results are integrated to produce the final classification results. This enables the RF classifier to have a higher tolerance for noise, and also avoids overfitting possibilities [63,64]. Moreover, RF has proven its capability of handling high-dimensional data by resulting in promising maps [31].

In this study, a pixel-based RF classifier was implemented within GEE because of the high potential of the RF classifier for mangrove ecosystem mapping [15,44,61]. In this regard, the final image mosaics, comprising seasonal Sentinel-1 and Sentinel-2 data, were inserted into an RF classifier. Meanwhile, half of the reference samples were employed to train the RF classifier. The RF classifier has several tuning parameters that affect the training phase of the classification step and, thus, directly influence the classification results. Among these parameters, the number of trees and variables at each node are the most influential parameters [63], which were set at 100 and equal to the square root of the number of input features, respectively. It should be noted that these values were determined through several trial and error attempts based on computational efficacy, visual inspection of the classification results, and the average out-of-bag sampling error of the RF classifier in the training phase.

#### *3.3. Accuracy Assessment*

Any thematic map derived from remote sensing data should be subjected to a trustworthy accuracy assessment to ensure its quality and reliability [65]. Therefore, both visual interpretation and statistical approaches were conducted to evaluate the accuracy of the final mangrove ecosystem map. High-resolution satellite images available within Google Earth and ArcMap were employed for visual assessment of the classification results. In terms of statistical accuracy assessment, independent test samples were incorporated to generate the confusion matrix, and then other metrics, such as OA, KC, producer accuracy (PA), and user accuracy (UA), were derived. The PA represents how well a specific area can be mapped, whereas the UA is an indicator of how well the produced map represents what really exists in the study area [66]. Furthermore, the F-score, the harmonic average of recall and precision, was computed for each class, and then its macro-averaging was reported for the classification results [67]. As already mentioned in Section 2.2.1, the classification procedure was repeated ten times with different training and test sets to comprehensively assess the classification performance for mangrove ecosystem mapping through a cross-validation step.

#### **4. Results**

Figure 3 presents the resulting mangrove ecosystem map (i.e., only classes within the mangrove ecosystem) with 10 m spatial resolution using RF classifier and the integration of seasonal SAR and optical satellite data within GEE. It should be noted that this map was produced based on a majority voting step of ten classification results. The majority voting step allows the production of a more reliable and accurate mangrove ecosystem map [68]. This is rooted in the fact that combining decisions of several classification results can lead to better recognition results [69]. Visually, the thematic map had an acceptable accuracy, indicating the high potential of the proposed method for delineating different classes. Generally, mangrove areas were depicted precisely, and their surroundings were also classified as aerial roots, as these roots grow around mangroves. Furthermore, the middle parts of the tidal channels were correctly classified as deep water, while other parts located near the coastline areas (i.e., tidal zone) were successfully distinguished as shallow water. Additionally, along with their corresponding high-resolution satellite images, two zoomed-in areas are also provided in Figure 3 for better visual interpretation. For instance, Figure 3b,c illustrates a zoomed-in area of the mangrove ecosystem, in which the mangrove areas are successfully identified. Moreover, based on Figure 3d,e, the proposed workflow was highly capable of discriminating between different classes with acceptable precision. In particular, mangrove areas and their surroundings (i.e., aerial roots) were properly discriminated, and also shallow water, deep water, and tidal zone classes were delineated appropriately. Furthermore, other areas without specific covers, which are wet (due to the existence of narrow water channels), were also accurately classified as mudflat.

The consistency images of six classes within the mangrove ecosystem were also computed based on the classification results of ten iterations to determine the robustness of the proposed classification procedure and the consistency of each pixel. It should be noted that the consistency images were produced based on the number of label assignments of each pixel to different classes. In other words, the consistency images present how many times a pixel was assigned to a specific class. The consistency values ranged between one and ten, in which one means very low consistency, while ten means highly consistent. Figure 4 shows the consistency images of six classes within the mangrove ecosystem. It was observed that in almost all cases, the classifier was able to assign consistent labels to each pixel in ten iterations. In particular, mangrove pixels were consistently labeled as mangrove with a high consistency rate of over 85% (i.e., considering pixels with values of eight to ten). Additionally, the tidal zone, deep water, and mudflat classes achieved reasonable rates of stability of over 70%. In contrast, the classes of shallow water and aerial roots obtained lower stability rates. In general, most inconsistencies occurred at boundaries and fragmented locations, where mixed pixels exist, suggesting the necessity of satellite images with a higher spatial resolution for more accurate mapping. For instance, the aerial roots had the highest rate of inconsistency of about 13% (i.e., considering pixels with values of one). This is probably rooted in the fact that these roots grow around mangrove areas and in mudflat conditions, and thus, the classifier may encounter difficulties in distinguishing them. This was more serious in distinguishing aerial roots from mudflat. Furthermore, another source of inconsistency was found between the shallow water and deep water classes, which is in fact associated with their greater similarity. The inconsistency revealed the effect of different training sets and also implied that incorporating the majority voting step could increase the reliability of the final thematic map.

**Figure 3.** Mangrove ecosystem map of the study area produced using the random forest (RF) algorithm and a combination of the Sentinel-1 and Sentinel-2 satellite images within the Google Earth Engine (GEE) platform (**a**), two zoomed areas of the produced map (**b**,**d**), and their corresponding high-resolution satellite images for better visualization (**c**,**e**).

Figure 5 presents the average values of PA and UA along with their standard deviation computed for ten iterations of producing mangrove ecosystem maps. The proposed classification framework obtained high average OA, KC, and F-scores of 93.23% (±1.1), 0.92 (±0.012), and 0.92 (±0.011), respectively. Furthermore, almost all classes achieved high PAs and UAs of over 90%, demonstrating the high potential and applicability of the proposed method for detailed mangrove ecosystem mapping. It is evident that the mangrove classes obtained significant average PA and UA of about 94.4% and 94.5%, indicating the applicability of the implemented approach for accurate mangrove delineation. In fact, the results proved that the proposed approach can not only be used for detailed mangrove ecosystem mapping but also can be effectively employed for accurate mangrove extent mapping. The highest and lowest accuracies were related to deep water and aerial roots, respectively. The aerial roots acquired moderate average PA and UA of approximately 77.6% and 78.2%, respectively, illustrating the challenging task of delineating aerial roots with high accuracy, which is also in accordance with [25]. Furthermore, Figure 5 demonstrates that the proposed approach obtained stable PAs and UAs in all ten iterations, since their standard deviation values ranged between 1.6% and 3.8% for all classes, except for the aerial roots, for which the standard deviation of PAs and UAs were 7.4% and 5.9%, respectively.

**Figure 4.** Consistency images of six classes (**a**) mangrove, (**b**) tidal zone, (**c**) deep water, (**d**) shallow water, (**e**) mudflat, and (**f**) aerial roots in the mangrove ecosystem. The consistency values ranged between one and ten, indicating very low and very high consistency, respectively.

**Figure 5.** Average producer accuracy (PA) and user accuracy (UA) of different classes along with their standard deviation (black error bars) computed from ten iterations.

The confusion matrices of the classifications were also investigated, to provide an in-depth statistical assessment in each iteration (see Figure 6). Generally, the proposed approach obtained acceptable accuracy since less confusion happened between classes, and confusion matrices were almost diagonal. However, there are several spots where some confusions are observable in Figure 6, such as the confusion between deep water and shallow water, which occurred in all iterations. Additionally, relatively high confusions were also found between aerial roots and other classes of mangrove and mudflat, which were the reasons for the moderate accuracy of aerial roots.

**Figure 6.** Confusion matrices of mangrove ecosystem maps for ten iterations.

An average confusion matrix (i.e., by averaging and rounding of the corresponding elements in 10 confusion matrices) is also provided in Table 3 to provide a comprehensive overview of all ten confusion matrices at a glance. Table 3 shows that the proposed approach obtained high average OA and KC values of 93.23% and 0.92, respectively, indicating the high potential of this approach for detailed mangrove ecosystem mapping. Moreover, the average omission error (OE) and commission error (CE) values were respectively 7.56% and 7.18%. The dominant confusion was associated with the aerial roots class, causing the highest OE and CE values of 23.06% and 21.92%, respectively, which are also in accordance with [25]. In particular, the two highest confusions occurred interchangeably between aerial roots/mudflat and aerial roots/mangrove. This was because these roots are grown in mudflat areas, exactly around mangroves. Therefore, it was challenging to separate aerial roots from two other classes due to the spectral/backscattering similarities and also the existence of mixed pixels at 10 m spatial resolution. Additionally, the second highest confusion was also observed between deep water and shallow water, which was actually because of their spectral/backscattering similarity.

**Table 3.** Average confusion matrix of mangrove ecosystem mapping using the random forest (RF) algorithm and a combination of the Sentinel-1 and Sentinel-2 satellite images processed within the Google Earth Engine (GEE) platform.


#### **5. Discussion**

#### *5.1. General Findings*

Mangrove ecosystems are routinely inundated and are located in inter-tidal zones [21]. Therefore, it is difficult to carry out field surveys to collect reference samples with global positioning system (GPS) devices [15]. The availability of highly accurate reference samples through field surveys is critical for reliable mapping; however, precise visual interpretation of high-resolution satellite images can compensate for this limitation. The proposed method was able to achieve high average OA and KC values, suggesting the suitability of applying this method when conducting fieldwork is difficult. It should be noted that generating reference samples even from high-resolution satellite imagery through precise visual interpretations can cause uncertainties and does not obviate the importance of in situ reference sample collection for more reliable mapping.

Despite the importance of mapping the mangrove extents, which have been conducted more frequently using either traditional or cloud computing-based approaches, fewer studies were devoted to producing detailed mangrove ecosystem maps. However, mapping relevant classes within a mangrove ecosystem (i.e., aerial roots and mudflat) is required to enhance condition assessment and also to suggest appropriate solutions for protection and conservation [25]. In particular, mapping the distribution and condition of aerial roots

are highly required for mangrove status monitoring since these roots provide different morphological and physiological adaptions to adverse environmental conditions [70]. For instance, these roots facilitate gas exchange to adapt mangroves to live in oxygen-poor locations [71], protect mangroves from flood effects and hurricanes [70], and mitigate shoreline erosion to prevent mangrove loss [72]. Furthermore, these roots support mangrove growth and accretion by stabilizing mudflat conditions [71,73]. Therefore, obtaining precise information about aerial roots helps to conduct more profound monitoring of mangroves and also to determine the locations for efficient reforestation programs.

In this study, the GEE cloud computing platform was employed for mangrove ecosystem mapping. GEE provides an unprecedented opportunity to employ dense time-series data from free-of-charge satellite images, and also it contains many built-in machine learning and image processing algorithms for satellite data manipulations [29,30]. Furthermore, the existence of ready-to-use data, along with other characteristics, contributes towards developing efficient methods with higher rates of automation and, thus, decreases the dedicated time for data acquisition, calibrating, and preprocessing [26]. Despite these advantages, GEE also has several limitations, such as restrictions on using numerous features and training samples (i.e., encountering limit exceed), which were also addressed in previous studies [32].

#### *5.2. Comparison with the Latest Global Mangrove Maps*

The obtained results indicated the high applicability of the proposed method for mangrove ecosystem mapping and also mangrove extent delineation. Therefore, the mangrove extent derived from the proposed method was compared with the recent available global mangrove extent map [74]. This global mangrove map was produced at 25 m (0.8 arc-second) spatial resolution using multi-source remote sensing data of the Advanced Land Observation Satellite (ALOS) Phased-Array L-band SAR (PALSAR) and optical data of Landsat-5 and Landsat-7 through an Extremely Randomized Trees (ERT) classifier [74]. Figure 7 presents the mangrove extent resultants of two products, along with the corresponding high-resolution satellite images for better visualization. Based on Figure 7, it is obvious that the proposed approach outperformed the global mangrove product. Clearly, based on Figure 7a–c, the proposed method was resistant to misclassification, probably of mudflat and aerial roots (i.e., generally unvegetated areas between mangrove patches). Likewise, Figure 7d–f also illustrates that the proposed method was successful at delineating narrow mangrove patches which occurred along the tidal zone areas. Despite the subtle differences in the classification algorithms and the considered classes in each workflow, the main discrepancies between the two maps were almost associated with the spatial resolution of remote sensing satellite data, which suggests the use of Sentinel-1 and Sentinel-2 data. Therefore, the synergistic use of Sentinel-1 and Sentinel-2 can definitely enhance the recent global mangrove extent mapping [53]. Although the utility of 10 m spatial resolution indeed enhanced the mangrove extent mapping, as the results suggest, it is still necessary to use satellite images with a higher spatial resolution to permit the production of more precise detailed mangrove ecosystem maps. This is rooted in the fact that higher spatial resolutions are necessary to discriminate between different classes that are located near each other, and to reduce the mixed pixel effect, such as separating aerial roots from two other classes of mangrove and mudflat. However, the current high cost of acquiring commercial high-resolution images is the main limiting factor.

**Figure 7.** High-resolution satellite images of two mangrove areas (**a**,**d**), mangrove extents from the recent global mangrove map (**b**,**e**) [74], and mangrove extents derived from the proposed method (**c**,**f**).

#### *5.3. Contribution of Multi-Source Remote Sensing Data*

As mentioned earlier, the Sentinel-1 SAR and Sentinel-2 optical data were combined to produce an accurate mangrove ecosystem map. It has been well acknowledged that a combination of optical and SAR remote sensing data could enhance the classification results of land cover mapping [75]. However, few studies have combined optical and SAR data for detailed mangrove ecosystem mapping, and they were mainly conducted based on optical data [25,76]. Therefore, here, the contribution of multi-source remote sensing data for detailed mangrove ecosystem mapping was evaluated based on two metrics of OA and KC. In this regard, the explained workflow (Section 3) was repeated with two other scenarios of using only optical or SAR data. The obtained results for the multi-source data processing were comprehensively discussed in the previous section, and the average OA and KC were 93.23% and 0.92, respectively. The average OA and KC achieved for the scenario, in which only optical data were incorporated, were 90.67% and 0.89, respectively. Although the accuracies based on the optical data were still satisfactory, the decline in the accuracies indicates the contribution of multi-source remote sensing data for more accurate mangrove ecosystem mapping. In particular, the use of only optical data resulted in over 1.5% (i.e., averaged over ten iterations) loss of class accuracies of each tidal zone, deep water, shallow water, and mudflat classes. Additionally, the results using only SAR data revealed a weak performance (i.e., OA = 62.98% and KC = 0.58); however, the SAR data was successful at delineating mangrove areas with a high accuracy of over 90% in almost all iterations. When using only SAR data, the lower accuracy was mainly associated with the low capability of separating different classes (i.e., except mangrove) within the mangrove ecosystem, with an average loss of 31.94% in class accuracies. The obtained results suggest that the inclusion of SAR data in the classification task, along with optical data, leads to producing a precise and accurate mangrove ecosystem map. This is mainly due to the fact that these two data sources can provide complementary information about the spectral and physical properties of different classes [77,78].

#### *5.4. Comparison with Annual Downscaling*

The seasonal downscaling approach was implemented for mangrove ecosystem mapping. In fact, the dense time-series SAR and optical data, acquired in 2019, were downscaled to seasonal features for further analysis. This approach takes the benefits of multi-temporal image analysis, as well as reducing the computational complexity of considering all available images without downscaling. Furthermore, since the mangrove ecosystems are routinely located in the inter-tidal zone and are inundated [21], the classifier may be influenced by the inundation intensity [79], biasing the mangrove extent and area. One solution is to use a single date image, in which the lowest tidal condition occurred; however, finding an image with such criteria may be challenging for several reasons, such as the presence of cloud and low temporal resolution of satellites. Therefore, it is recommended to apply time-series data to acquire more reliable classification results. In this study, seasonal features were used; however, other researchers applied annual downscaling for mangrove mapping [39]. Consequently, the usefulness of seasonal downscaling is compared with the annual downscaling approach. In this regard, the mentioned workflow was also applied to annual downscaled SAR and optical data, and the obtained results were compared with the seasonal approach. As is clear from Figure 8, incorporating seasonal downscaling achieved higher OAs and KCs in all iterations and, thus, it suggests that the seasonal downscaling can provide more discriminative information compared to annual downscaling [80]. This may be rooted in the fact that annual downscaling would override the tidal fluctuations, while the seasonal downscaling includes seasonal variations of water level (i.e., dry and wet seasons), reflecting tidal fluctuations [78].

**Figure 8.** Overall accuracy (OA; (**a**)) and Kappa coefficient (KC; (**b**)) of annual and seasonal downscaling of SAR and optical data for mangrove ecosystem mapping.

#### **6. Conclusions**

Producing detailed mangrove ecosystem maps is essential for natural resource monitoring and sustainable development goals tracking. In this paper, a straightforward yet robust workflow was proposed to produce mangrove ecosystem maps with high accuracies. For this purpose, Sentinel-1 and Sentinel-2 data were inserted into an RF classifier within the GEE cloud computing platform. The final classification results were comprehensively evaluated through different analyses, and further discussions were provided. The classification results obtained high average OA and KC of 93.23% and 0.92, respectively. Moreover, all the classes obtained high accuracies, except the aerial roots, which achieved moderate accuracies in all cases, suggesting the consideration of other possible solutions (i.e., higher resolution images) for more accurate delineation of this class. Altogether, taking both qualitative (i.e., visual interpretation) and quantitative (i.e., statistical accuracy assessment) evaluation criteria into account, the proposed method confirmed its applicability in producing a detailed mangrove ecosystem map. Furthermore, the comparison results proved the contribution of multi-source remote sensing data (i.e., SAR + optical), as well as the effectiveness of seasonal downscaling. The proposed workflow was implemented based on an open-source platform and free-of-charge satellite data, making it appealing and applicable in almost all countries.

**Author Contributions:** Conceptualization, A.G. and S.Z.; methodology, A.G. and S.Z.; software, A.G.; validation, A.G., S.Z. and R.M.A.; formal analysis, A.G.; investigation, A.G.; data curation, A.G. and S.Z.; writing—original draft preparation, A.G. and R.M.A.; writing—review and editing, A.G., S.Z., R.M.A., M.A., S.J., A.M.; visualization, A.G.; supervision, A.G., A.M. and S.J.; project administration, A.G.; funding acquisition, M.A. and S.J. All authors have read and agreed to the published version of the manuscript.

**Funding:** S.J. was partially funded by Knut and Alice Wallenberg grant 99007-RFh2018-2016.

**Data Availability Statement:** The produced mangrove ecosystem map using Sentinel-1, Sentinel-2, and Random Forest algorithm and the derived mangrove extent map are available upon request from the first author.

**Acknowledgments:** The authors would like to thank the European Space Agency (ESA) to freely provide Sentinel-1 and Sentinel-2 satellite data. Additionally, the authors thank the anonymous three reviewers for their comments.

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

#### **References**


## *Article* **Global Sensitivity Analysis for Canopy Reflectance and Vegetation Indices of Mangroves**

**Chunyue Niu \*, Stuart Phinn and Chris Roelfsema**

Remote Sensing Research Centre, School of Earth and Environmental Sciences, The University of Queensland, Brisbane, QLD 4072, Australia; s.phinn@uq.edu.au (S.P.); c.roelfsema@uq.edu.au (C.R.)

**\*** Correspondence: c.niu@uqconnect.edu.au

**Abstract:** Remote sensing has been applied to map the extent and biophysical properties of mangroves. However, the impact of several critical factors, such as the fractional cover and leaf-to-total area ratio of mangroves, on their canopy reflectance have rarely been reported. In this study, a systematic global sensitivity analysis was performed for mangroves based on a one-dimensional canopy reflectance model. Different scenarios such as sparse or dense canopies were set up to evaluate the impact of various biophysical and environmental factors, together with their ranges and probability distributions, on simulated canopy reflectance spectra and selected Sentinel-2A vegetation indices of mangroves. A variance-based method and a density-based method were adopted to compare the computed sensitivity indices. Our results showed that the fractional cover and leaf-to-total area ratio of mangrove crowns were among the most influential factors for all examined scenarios. As for other factors, plant area index and water depth were influential for sparse canopies while leaf biochemical properties and inclination angles were more influential for dense canopies. Therefore, these influential factors may need attention when mapping the biophysical properties of mangroves such as leaf area index. Moreover, a tailored sensitivity analysis is recommended for a specific mapping application as the computed sensitivity indices may be different if a specific input configuration and sensitivity analysis method are adopted.

**Keywords:** global sensitivity analysis; PAWN; canopy reflectance model; vegetation index (VI); mangroves

#### **1. Introduction**

The extent and biophysical properties of mangroves have been mapped from multispectral remote sensing (RS) images [1–4]. However, several factors such as the woody material and tidal height could affect the canopy reflectance [5,6], which was often ignored. The scientific community has been aware of these, but quantitative assessments are still needed to identify essential factors that have a major impact on the spectral responses of mangroves.

Sensitivity analysis (SA) is the study of how to attribute the variation or uncertainty in the output of a model to variations in the model inputs [7,8]. It can be employed to rank the input factors based on their contributions to the variations of the model output, to identify factors that have an insignificant impact on output variations and to find out regions in the input space that produce particular output values [7,8]. SA can be divided into global SA (GSA) or local SA according to whether input values are allowed to vary across the entire input space or around a reference point in the input space [8].

SA has been broadly employed in RS to recognise the input factors that had a major impact on vegetation canopy reflectance using various canopy reflectance models (CRMs). The PROSPECT model [9], the SAIL family including SAILH, 4SAIL2 and soil–leaf–canopy (SLC) [10–12] and the PROSAIL (PROSPECT + SAIL) model [13] were among the CRMs that had been extensively explored. Variance-based SA (VBSA) methods have attracted a great deal of attention. Its fundamental concept is to decompose the direct contribution of a factor to output variance and its overall contribution when interaction effects with

**Citation:** Niu, C.; Phinn, S.; Roelfsema, C. Global Sensitivity Analysis for Canopy Reflectance and Vegetation Indices of Mangroves. *Remote Sens.* **2021**, *13*, 2617. https:// doi.org/10.3390/rs13132617

Academic Editors: Lars T. Waser and Chandra Giri

Received: 13 May 2021 Accepted: 20 June 2021 Published: 3 July 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

other factors are considered [8,14,15]. Corresponding metrics are usually termed as firstand total-order sensitivity indices. Liu et al. [16] applied the extended Fourier amplitude sensitivity test (EFAST) [17] to calculate the first order indices and interaction effects for analysing the sensitivities of vegetation indices (VIs) to the leaf area index (LAI) and other interference factors, such as leaf chlorophyll content (*C*ab) and average leaf inclination angle (*θ*l), and evaluating the uncertainties caused by these factors. This study found that leaf parameters such as *C*ab and leaf dry matter content (*C*m) had an increasing impact on the uncertainty of estimated LAI towards higher LAI and *θ*<sup>l</sup> was the most critical factor for LAI estimation if an ellipsoidal distribution [18] was used. Xiao et al. [19] also adopted EFAST and model simulations to investigate the sensitivities of reflectance and VIs to biophysical and biochemical parameters at the leaf, canopy and regional levels. The importance of leaf parameters reduced at the canopy level, especially for a LAI value of 0–3, where LAI dominated except for the absorption bands of *C*ab and leaf water content (*C*w). Although the contributions of LAI and soil were still evident for a sparse canopy, the fractional cover of vegetation (fCv) dominated at the regional level. The importance of LAI and soil dropped significantly for LAI > 3 at both the canopy and regional levels. Mousivand et al. [20] proposed an improved design and sampling for SA to identify influential and non-influential factors on canopy reflectance based on the SLC model [12], in which soil moisture (SM) and fCv were explicitly considered. fCv, LAI, *θ*<sup>l</sup> and SM were recognised as the most influential factors for a LAI range of 0–6.

The studies mentioned above mainly focused on terrestrial vegetation, while only limited SA studies were carried out for aquatic vegetation. Villa et al. [21] introduced two new VIs specifically for aquatic vegetation: the normalised difference aquatic vegetation index (NDAVI) and the water adjusted vegetation index (WASI) and performed VBSA to evaluate their sensitivities to LAI, *C*ab and *θ*l. Both radiative transfer simulations and linear spectral mixture simulations based on real endmembers were applied. The new VIs were found to have a higher sensitivity to LAI and *θ*<sup>l</sup> and better capacity in separating aquatic and terrestrial vegetation compared with previous VIs for terrestrial vegetation such as NDVI. Unfortunately, the fCv in 4SAIL2 was fixed to 1, and thus its impact was not discussed. Zhou et al. [22] conducted EFAST for emergent and submerged vegetation based on a CRM model for aquatic vegetation [23]. Four scenarios were adopted according to the combinations of shallow or deep water and sparse or dense canopies. Canopy reflectance spectra and four VIs, including NDAVI and WAVI for Sentinel-2A (S2A) images, were simulated. The most influential factors were found to be different for two vegetation types, with leaf and canopy parameters being dominant for emergent vegetation while water components accounted for most variability in canopy reflectance of submerged vegetation. For emergent vegetation, the total order sensitivity indices of all four VIs to LAI were high for a sparse canopy but dropped notably for a dense canopy. NDAVI outperformed the other three VIs for a dense emergent canopy in terms of their sensitivities to LAI.

Although extensively applied, VBSA was based on an implicit assumption that variance was sufficient to represent output variability [24]. This assumption was not appropriate if the output was highly skewed or had multiple modals [25]. Moreover, the probability distribution of an input factor also mattered, and the total order indices could not be computed if there were correlated inputs [24]. However, input probability distributions were usually assumed to be uniform, e.g., [20,22]; the potential correlations between input factors and the output distributions were rarely examined in many VBSA studies.

These limitations of VBSA could be avoided by adopting moment-independent sensitivity indices that relied on the entire output distribution rather than a particular moment [26]. These were also referred to as density-based SA (DBSA) as they employed the probability density function (PDF) or cumulative distribution function (CDF) of the output to characterise its uncertainty [8,27]. The basic idea was to evaluate the sensitivity by calculating the divergence between unconditional output PDF and conditional output PDF, where "unconditional" represented varying all input factors while "conditional" meant that the corresponding input factor was fixed [8]. Despite its advantages, DBSA was not

widely applied, partly due to the challenges in obtaining conditional PDFs. Thus, using CDF instead of PDF could reduce the computational costs and make it relatively easy to implement [25].

In summary, current studies on SA of canopy reflectance of terrestrial vegetation and their conclusions might not be applied to mangroves directly, partly because the adopted CRMs ignored essential factors for mangroves such as fCv, water body and woody material. Additionally, a few recently proposed VIs for aquatic vegetation, e.g., [21], would possibly be effective for the RS of mangroves, but quantitative assessments were required. In addition, the impact of potential correlations between input factors, e.g., plant area index (PAI) and fCv [28], input probability distributions and the adopted SA methods on the produced sensitivity indices also needed more evaluations. Therefore, the objective of this study is to perform a systematic GSA for mangroves to quantitatively assess the impact of various biophysical and environmental factors on their canopy reflectance and selected VIs, and to identify the most influential factors under different scenarios, e.g., sparse or dense canopies. The findings may help to improve our understanding of the spectral characteristics of mangroves and recognise the factors that may be mapped from multispectral satellite images.

#### **2. Methodology**

#### *2.1. Overview*

In this study, 16 factors were selected based on a priori knowledge, similar studies such as [19,20,22] and trial tests. As only images captured by nadir or near-nadir viewing satellites, e.g., Sentinel-2, would be applied at this stage, the observation zenith angle was set as 0◦, and the solar zenith angle was fixed as 30◦ to avoid distraction. Then various mangrove scenarios such as sparse or dense canopies were set up according to one or more properties such as the crown PAI and fCv. This corresponded to different mangrove formations such as open or closed forests. Samples of the 16 input factors were generated using VBSA or DBSA and following uniform or normal probability distributions. The canopy reflectance spectra from 400 to 2500 nm with a 1 nm interval were then simulated using a CRM of mangrove [28]. Additionally, values of corresponding S2A bands (Bands 2–8, 8a, 11 and 12) were simulated by convolving the simulated reflectance spectra with the spectral response functions of S2A. Then twelve VIs were calculated. Finally, both VBSA and DBSA were performed and the computed sensitivity indices were examined and compared. The GSA flowchart is shown in Figure 1.

#### *2.2. Canopy Reflectance Model of Mangroves*

The adopted CRM of mangroves was a one-dimensional multiple-layer radiative transfer model [28]. It considered essential biophysical and biochemical properties of mangroves such as *C*ab, PAI, fCv, L2T ratio, inclination distributions of leaves and woody material and environmental factors such as the water body. It could be regarded as an extension of a previous CRM for aquatic vegetation [23], which was a revision of the PROSAIL model. Mangroves were divided into three layers as in [28]: the understory (L1), the stem (L2) and the crown (L3) layers (Figure 2). The layer index would be added to the variable names in the following to distinguish them from each other, e.g., PAI(3) for the PAI of the crown layer. The height of the understory was assumed to be the same as water depth (*H*w), and water was not allowed to immerse the crown as that could change canopy reflectance significantly [29]. Stems were assumed to be covered by the crowns and hence only had a minor contribution to canopy reflectance. Thus, the stem layer was assigned a low PAI value of 0.5 and a low fCv value of 0.15, according to [28]. Similarly, most factors with regard to the optically active components in the water body and other environmental factors such as bottom reflectance were also fixed as their impact on canopy reflectance was minor when covered by the canopy. For example, the contribution of suspended particles in the background water on canopy reflectance was around 10% or lower and mainly in the visible region [22,28]. Therefore, more attention was paid to the crown, the understory

and *H*w, and other factors were beyond the scope of this study. *H*<sup>w</sup> corresponded to the changing tidal height. PAI(1) and fCv(1) were employed to explore if there were any potential contributions from the understory to the observed canopy reflectance, especially for a sparse canopy.

**Figure 1.** The flowchart of global sensitivity analyses for mangroves.

**Figure 2.** Illustrations of the understory (L1), stem (L2) and crown (L3) layers for mangroves. *H*w is the water depth and is assumed to be the same as the height of the understory for simplicity (The symbols are courtesy of the Integration and Application Network, University of Maryland Center for Environmental Science (ian.umces.edu/symbols/, accessed on 14 June 2018)).

#### *2.3. Mangrove Scenarios*

If not specified, the input factors were independent, as they were treated in other SA studies [19,20,22]. All samples discussed in Section 2.3.1 were generated using the Latin hypercube sampling (LHS) [30] following a uniform or normal probability distribution. Although independent input samples were widely adopted in SA studies, there might be an issue with some biophysical properties. For example, it was possible to obtain a high PAI(3) value together with a low fCv(3) value or vice versa if they were assumed to be independent (Figure 3a). However, this was not reasonable in nature and the field data

in [28] revealed a positive correlation with a coefficient of determination of 0.73 between the PAI and fCv of mangroves. Therefore, another input dataset with correlated PAI(3) and fCv(3) was also designed in Section 2.3.2.

**Figure 3.** Illustrations of generated samples (*N* = 1000) if the crown plant area index (PAI(3)) and fractional cover (fCv(3)) are assumed to be independent (using Latin hypercube sampling) (**a**) or linearly correlated (**b**).

#### 2.3.1. Factors, Their Ranges and Distributions

GSA was performed for three different mangrove scenarios defined by the ranges of PAI(3) and fCv(3): general canopy (GEN), sparse canopy (SPS) and dense canopy (DEN) (Table 1), as previous studies demonstrated that different ranges of inputs could result in different sensitivity values [19,22]. PAI(3) ranges were set as 0–6 for GEN, 0–3 for SPS and 3–6 for DEN according to the literature [19,31]. Correspondingly, fCv(3) ranges were 0–1, 0–0.7 and 0.5–1. The overlap between the fCv(3) ranges for SPS and DEN was to reflect the variability in mangrove forests. The SPS and DEN scenarios corresponded to open and closed mangrove forests, while the GEN scenario included both situations and implied a significant variation of mangrove formations. Ranges of other factors such as leaf parameters were not changed across these scenarios.



Note: Uniform or normal is the probability distribution of input factors, with Min and Max for minimum and maximum values and Mean and STD for mean value and standard deviation. <sup>a</sup> Equivalent to average inclination angle (AIA) by a = (45 − AIA)∗π2/360 [11].

In addition to the GSA methods and ranges of inputs, the probability distributions of inputs might also impact the sensitivity values [32]. The uniform (UNI) and normal (NRM) probability distributions were both adopted for the input factors under SPS and DEN. The mean values for NRM were selected according to [28], and corresponding standard deviation values were set as 10% of the mean values, except for LIDFa(3) and WIDFa(3). NRM indicated that a priori information on the mangrove forest was available and the variability in canopy structure was small, while UNI indicated a larger variation in canopy structure. Under GEN, only a uniform distribution was applied as it might not be suitable to use a nominal value to represent such a wide range of PAI(3).

#### 2.3.2. Correlated PAI(3) and fCv(3)

In this case, the linear relationship between PAI and fCv of mangroves from [28] was applied to generate samples for correlated PAI(3) and fCv(3). First, the ranges of PAI(3) and fCv(3) were set as 1–5 and 0.1–0.9. Hence, their values were within 0–6 and 0–1 after adding random variations. Then random samples of all 16 factors were generated as in the GEN case. Second, PAI(3) or fCv(3) was modified to make them correlated. For VBSA, fCv(3) values were recalculated according to PAI(3) values and the linear relationship, and subsequently, random variations within ±0.1 were added (Figure 3b). In the case of PAWN, the same methodology was used to compute unconditional CDF and conditional CDFs for most factors except fCv(3). When fCv(3) was fixed for computing the corresponding conditional CDFs, PAI(3) was recalculated based on fCv(3), and then a random variation within ±0.8 was given. Finally, sensitivity indices for these two sample datasets were computed using VBSA and PAWN, respectively. Correlated inputs were only applied to GEN as a demonstration because the relationship between PAI and fCv could be different for various sites and species.

#### *2.4. Vegetation Indices*

VIs have been broadly applied for mapping the biophysical properties of mangroves such as LAI [2,3,31], but the impact of woody material, fCv and background water was rarely reported. Although a large number of VIs were available [33] (https://www.indexdatabase.de/, accessed on 9 December 2019), this study did not intend to investigate all possible choices. Instead, the VIs in Table 2 were simulated based on the spectral response function of the S2A multispectral instrument to assess their sensitivities to the input factors such as PAI, fCv, L2T ratio and *H*w. These indices were either widely applied for terrestrial vegetation such as NDVI [34] and EVI [35] or in particular, proposed for aquatic vegetation such as NDAVI [36]. Some of the VIs for aquatic vegetation included shortwave infrared (SWIR) bands to address the existence of background water, e.g., RGVI [37] and WFI [38].


**Table 2.** Definitions and sources of the adopted vegetation indices in this study.

#### *2.5. Sensitivity Analysis Methods*

#### 2.5.1. Variance-Based Sensitivity Analysis

VBSA employed variance to represent sensitivity and attributed the variance of the model output to the variances of input factors. For a model with *k* independent input factors and a scalar output *Y*: *Y* = *f*(*X*1, *X*2, ... , *Xk*), the variance (*V*(·)) of *Y* can be decomposed as [14,15]:

$$V(Y) = \sum\_{i} V\_{i} + \sum\_{i} \sum\_{j>i} V\_{i\bar{j}} + \dots + V\_{12\dots k} \tag{1}$$

The first order indices *S*<sup>i</sup> could be represented by the expected reduction in output variance if *Xi* could be fixed [15]:

$$S\_i = \frac{V\_i}{V(Y)} = \frac{V\_{X\_i}(E\_{X\_{\sim i}}(Y|X\_i))}{V(Y)} \tag{2}$$

where *E*(·) was the expectation operator and subscripts *Xi* or **X**~*<sup>i</sup>* indicated that the operation (*E*(·) or *<sup>V</sup>*(·)) was taken over *Xi* or all factors but *Xi*. The total order indices *<sup>S</sup>*<sup>T</sup> *<sup>i</sup>* could be written as [15]:

$$S\_i^T = \frac{E\_{\mathbf{X}\_{\sim i}}(V\_{\mathbf{X}\_i}(Y|\mathbf{X}\_{\sim i}))}{V(Y)} = 1 - \frac{V\_{\mathbf{X}\_{\sim i}}(E\_{\mathbf{X}\_i}(Y|\mathbf{X}\_{\sim i}))}{V(Y)} \tag{3}$$

Likewise, *V***X**∼*<sup>i</sup>* (*EXi* (*Y* **X**∼*i*)) was the expected reduction in output variance if all input factors but *Xi* could be fixed. Other partial variances in Equation (1) corresponded to the interactions among input factors, and more details could be found in [15]. Both first and total order indices were broadly used, but this study only focused on total order indices. Due to numerical errors, this approach might produce negative sensitivity indices, usually corresponding to noninfluential factors [7]. Hence, negative values were reset to zero.

Open-source software was available for VBSA, such as SimLab [43] and the SAFE toolbox [44] (https://www.safetoolbox.info/info-and-documentation/, accessed on 7 May 2019). The SAFE toolbox for MATLAB® was employed in this study as it also provided a DBSA method. Trial tests showed that the computed sensitivity indices in VBSA could be affected by the number of samples, especially for NRM. Multiple sample sizes from 16,000 (1000 × 16) to 128,000 (8000 × 16) with an increment of 16,000 were tested for SPS\_NRM and DEN\_NRM. Finally, 112,000 samples were generated for each VBSA case as the computed sensitivity indices of reflectance spectra were consistent and stable.

#### 2.5.2. Density-Based Sensitivity Analysis

The PAWN method [25] in the SAFE toolbox was adopted for DBSA in this study. It employed output CDFs instead of PDFs to reduce computational costs and to avoid tuning parameters, e.g., the bin width for calculating PDF [25]. The distance between unconditional and conditional output CDFs was taken as a measure of sensitivity. First, *N*<sup>u</sup> samples were randomly generated for a model of *k* input factors *Y* = *f*(*X*1, *X*2, ... *Xk*), to evaluate the empirical unconditional CDF of output *FY*(*Y*). Second, *n* samples were generated for each input factor *Xi* and acted as conditional points. At each conditional point, *N*<sup>c</sup> random samples were generated over **X**~*<sup>i</sup>* to obtain the empirical conditional CDF of output *FY*|*Xi* (*Y*). The Kolmogorov–Smirnov (KS) statistic was applied as a measure of the distance between unconditional and conditional CDFs (as cited in [25]):

$$\text{KS}(X\_{\bar{i}}) = \max\_{\mathcal{Y}} \left| F\_{\mathcal{Y}}(\mathcal{Y}) - F\_{\mathcal{Y}|X\_{\bar{i}}}(\mathcal{Y}) \right| \tag{4}$$

In other words, it was the maximum divergence between conditional and unconditional CDFs for a conditional point of *Xi*. Considering all the conditional points of an input factor *Xi*, another statistic, e.g., the maximum or the median, was computed over all conditional points of *Xi* to obtain the PAWN sensitivity index *P*<sup>T</sup> *<sup>i</sup>* [25]:

$$P\_i^T = \text{statistic}[KS(X\_i)]\tag{5}$$

By its definition, *P*<sup>T</sup> *<sup>i</sup>* would be between 0 and 1, and higher values suggested a greater influence on the output [25]. Taking a result for S2A NDVI as an example, the unconditional and conditional CDFs and the KS statistics for all conditioning points of fCv(3) are in Figure 4. As the KS curves had various patterns, derived *P*<sup>T</sup> *<sup>i</sup>* values using the maximum or the median statistic could be different and thus lead to divergent conclusions. The maximum statistic was adopted because it could detect factors that had influences on the output for at least one conditioning point [32]. The number of samples in a PAWN case was 66,000 in this study, with 2000 samples used to calculate the unconditional CDF, and 20 conditioning points for each of the 16 input factors and 200 samples used per conditioning point. More conditioning points or samples did not change the computed sensitivity indices remarkably.

**Figure 4.** Unconditional (red) and conditional (grey) cumulative distribution functions (**a**) and corresponding Kolmogorov– Smirnov (KS) statistics (**b**) of the normalised difference vegetation index with regard to the fractional cover of the mangrove crown (fCv(3)). The PAWN index can be different if the maximum (blue) or the median (green) metric is adopted (**b**).

#### **3. Results**

The normalised sensitivity indices (NSIs) for reflectance spectra and sensitivity indices (SIs) for S2A bands and VIs from VBSA and PAWN were presented as stacked bar plots (Figures 5 and 6). As VIs had various scales, they were not normalised to avoid confusion, especially for VBSA indices that sometimes only one or two factors achieved positive sensitivity values, e.g., MDI2 in Figure 6c.

**Figure 5.** Stacked bar plots for normalised sensitivity indices (NSIs) of reflectance spectra in VBSA (**a**–**e**) and PAWN (**f**–**j**) for a general (GEN) canopy with uniformly sampled inputs (UNI) (**a**,**f**), a sparse (SPS) canopy with UNI (**b**,**g**), SPS with normally sampled inputs (NRM) (**c**,**h**), a dense (DEN) canopy with UNI (**d**,**i**) and DEN with NRM (**e**,**j**). Each *x* value corresponds to a wavelength. The NSI value is represented by the height of a bar and marked in the *y*-axis. Each colour corresponds to a factor defined in Table 1.

**Figure 6.** Stacked bar plots for sensitivity indices (SIs) of Sentinel-2A (S2A) band reflectance and vegetation indices (VIs) in VBSA (**a**–**e**) and PAWN (**f**–**j**) for a general (GEN) canopy with uniformly sampled inputs (UNI) (**a**,**f**), a sparse (SPS) canopy with UNI (**b**,**g**), SPS with normally sampled inputs (NRM) (**c**,**h**), a dense (DEN) canopy with UNI (**d**,**i**) and DEN with NRM (**e**,**j**). Each *x* value corresponds to a S2A band reflectance or a VI. The SI value is represented by the height of a bar and marked in the *y*-axis. The scale of the *y*-axis is restricted from 0 to 3.8 to show sufficient details. Each colour corresponds to a factor defined in Table 1.

#### *3.1. General Scenario*

For GEN, two large areas in blue and yellow in Figure 5a could be easily recognised, indicating that L2T(3) and fCv(3) had significant influences on canopy reflectance across nearly the whole visible to SWIR regions in VBSA. Exceptions were the boundaries of chlorophyll absorption bands around 600 nm and 710 nm where the impact of *C*ab rose and fell abruptly and the NIR bands where the influence of L2T(3) dropped. WIDFa(3) was the next most influential factor, especially in the SWIR, followed by PAI(3) across NIR and SWIR bands. Comparatively, PAWN NSIs were slightly different from VBSA NSIs in that PAI(3) had higher NSI values while the NSI values of L2T(3) and fCv(3) were lower, especially in SWIR (Figure 5f). It is worth noting that *H*w had higher NSI values in visible bands compared with its NSI values in other regions.

Since PAWN SIs were always between 0 and 1 while VBSA SIs could have various scales and negative values were reset to zero, the NSI patterns of VIs in the two methods could appear distinct. The stacked PAWN SIs were around 3.5 while stacked VBSA SIs varied from lower than 1 to higher than 6. The scale of the *y*-axis in Figure 6 was limited to 0–3.8 to maintain sufficient details. Once they were normalised, the length of each bar could be stretched or compressed. Thus, the relative length of the bar for a specific factor might look different in stacked bar plots of SIs and NSIs. In PAWN, multiple VIs had similarly high SI values of fCv(3) such as S2A-B11, MNDWI and RGVI (Figure 6f). A challenge with the PAWN indices was that all factors obtained a SI value between 0 and 1. Hence, the influential factors might look less outstanding than they did in VBSA, even if its SI value was high, e.g., the NSI of fCv(3) for MNDWI was 0.21, but its SI value was 0.68.

The output PDFs in VBSA were shown in Figure 7. Those outputs that had high NSI or SI values to fCv(3) in VBSA, such as S2A-B8a, MNDWI, NDAVI and RGVI, exhibited various PDFs. Their SIs with confidence intervals were obtained via bootstrapping and further examined. Although the PDF of RGVI was remarkably skewed, all factors had narrow confidence intervals (around 0.1 or lower) and the most influential factor fCv(3) was outstanding (not shown). The confidence intervals were also narrow for MNDWI and NDAVI but were relatively large for S2A-B8a (around 0.25 or higher, Figure 8a). In the case of PAWN, the KS plots for MNDWI with regard to input factors in Figure 8b could provide more insights. KS values of fCv(3) would increase if fCv(3) was approaching 0 or 1 from 0.4 and for the latter case, KS values would become stable after fCv(3) reached about 0.9. When L2T(3) moved towards the two ends 0 or 1, its KS values slightly increased but were not obvious. As for PAI(3) and *H*w, their KS values only grew noticeably if their values decreased to zero, e.g., lower than 1 for PAI(3). This implied that PAI(3) was more influential when its value was low.

**Figure 7.** Probability density function (PDF) of Sentinel-2A (S2A) band reflectance and vegetation indices for a general (GEN) canopy with uniformly sampled inputs (UNI) (blue curves), a sparse (SPS) canopy with UNI (red curves) and SPS with normally sampled inputs (NRM) (yellow curves). The bin size for SPS\_NRM is one-third of that for GEN\_UNI and SPS\_UNI to show details of the latter two.

**Figure 8.** Mean sensitivity indices and confidence intervals estimated via bootstrapping for S2A-B8a in VBSA (**a**). Kolmogorov–Smirnov (KS) statistics for MNDWI with regard to input factors in PAWN (**b**).

#### *3.2. Sparse Mangroves—Uniform Input Probability Distributions*

The NSIs of reflectance spectra for SPS\_UNI had many similarities with those for GEN\_UNI, but there was an increased impact of *H*<sup>w</sup> in visible bands (Figure 5b,g). In the NIR region, the NSI values of PAI(1) and fCv(1) slightly increased in PAWN, but only the NSI values of fCv(1) increased in VBSA. This implied that the understory and water background had a larger impact on SPS, which was intuitive. Therefore, more attention was paid to *H*w in this section.

S2A-B3 achieved the highest SI value of 0.87 to *H*w, but all factors had large confidence intervals (around 0.4 or higher, Figure 9a). Although *H*<sup>w</sup> could still be distinguished, its bounds were already overlapped with those of other factors. This was somewhat unexpected as the S2A-B3 values were quite normally distributed and far from skewed (Figure 7c). Moreover, a larger sample size of up to 128,000 did not significantly reduce the confidence intervals. For RGVI and NDAVI, *H*w, fCv(3) and PAI(3) were all identified as influential, but their SI values slightly varied and their confidence intervals might overlap. For NDWI, the SI of fCv(3) was slightly higher than that of *H*w, i.e., 0.360 vs. 0.348. Their confidence intervals were narrow (lower than 0.1) but still overlapping. Ranking them further might not be reliable because of the overlapping confidence intervals [8,45]. Therefore, a high SI itself in VBSA might not guarantee a solid link between the output and an input factor. It was better to check the output distributions and confidence intervals as well, as suggested by [32].

**Figure 9.** Mean sensitivity indices and confidence intervals estimated via bootstrapping for S2A-B3 in VBSA (**a**). Kolmogorov–Smirnov (KS) statistics for NDWI with regard to input factors in PAWN (**b**).

By comparison, S2A-B3 had similar SI patterns in PAWN, with the SI of *H*<sup>w</sup> higher than 0.6 and those of other factors lower than 0.3, but the confidence intervals of all factors were much narrower (around 0.1 or lower). For NDWI, *H*w, fCv(3) and PAI(3) had SI values of 0.685, 0.414 and 0.404, separately. All factors had narrow confidence intervals (lower than 0.1) in PAWN, and there was no overlapping between them and other factors. As for the KS plot Figure 9b, PAI(3) and fCv(3) presented similar patterns with those in Figure 8b, except that the KS values of fCv(3) continued rising with increasing fCv(3) since its maximum value was 0.7 rather than 1 in GEN. *H*w had a two-way impact on NDWI and S2A-B3. This might be attributed to the different contributions of absorption and scattering by the water body and its components, both of which could vary with water depth. In contrast, only low *H*<sup>w</sup> values could make a difference for VIs based on SWIR bands, e.g., MNDWI in Figure 8b, as the absorption of water in SWIR was strong.

#### *3.3. Sparse Mangroves—Normal Input Probability Distributions*

This case corresponded to small variations in the biophysical and biochemical properties of mangrove plots and the environmental factors. As shown in Figure 5b,c, there was a dramatic change in the NSI patterns of reflectance spectra in VBSA. fCv(3) and LIDFa(3) became dominant in NIR, while the most influential factor was L2T(3) in SWIR. Additionally, the impact of other factors was minor in these regions. Consequently, only a limited number of factors obtained positive SIs for multiple VIs such as the EVI and NDVI (Figure 6c). In contrast, the NSI patterns in PAWN also changed but more smoothly (Figure 5h). fCv(3) and L2T(3) were also the most influential factors in NIR and SWIR, separately. In addition, the influences of fCv(1) and LIDFa(3) increased in NIR and the impact of tree shape factor zeta(3) also increased.

As shown in Figure 7, the output PDFs in VBSA for SPS\_NRM were remarkably different from those for GEN\_UNI and SPS\_UNI. This indicated that input distributions could affect the output distributions. However, the normal-like output PDFs did not guarantee that the confidence intervals of SIs were narrow. Taking MNDWI as an example, the confidence intervals were around 0.3, and there was overlapping between influential factors fCv(3) and L2T(3) and other factors (Figure 10a). Comparatively, the confidence intervals for MNDWI (Figure 10b), MDI1 and LSWI in PAWN were basically within 0.1.

**Figure 10.** Mean sensitivity indices and confidence intervals estimated via bootstrapping for MNDWI in VBSA (**a**) and PAWN (**b**).

#### *3.4. Dense Mangroves—Uniform Input Probability Distributions*

When the canopy became dense, the impact of PAI(3) dropped quickly, and the impact of fCv(3) also declined noticeably, as shown by their shrunken areas in Figure 5d,i. *H*<sup>w</sup> barely made a difference in this case. In the meanwhile, the influences of LIDFa(3) increased in NIR and *C*w in SWIR. Moreover, the impact of L2T(3) and WIDFa(3) grew significantly, especially in SWIR. L2T(3) had the highest SI values over other factors for multiple VIs such as MDI2, LSWI, NDAVI and WFI (Figure 6d,i). This might be attributed to the distinct spectral responses of leaves and woody material, and the adopted reflectance values of woody material were higher than the reflectance of leaves in SWIR.

It was interesting to find that *C*ab also had outstanding PAWN SI values for multiple VIs such as NDVI, EVI and S2A-B4 (Figure 6i). Comparatively, this was less significant for NDVI and EVI in VBSA (Figure 6d). The KS plots for NDVI in PAWN were presented in Figure 11b. L2T(3) had a two-way impact while the KS values of *C*ab increased dramatically only when its values were lower than 20 <sup>μ</sup>g·cm<sup>−</sup>2. This was because low *<sup>C</sup>*ab values close to zero could remarkably increase the canopy reflectance in chlorophyll absorption bands [46]. The SI values of *C*ab and L2T(3) for NDVI in PAWN were 0.879 and 0.771, respectively, with SI values of other factors lower than 0.4. The confidence intervals of all factors were within 0.1 in PAWN. In contrast, L2T(3) was the most influential factor (SI = 0.544) for NDVI in VBSA, followed by *C*ab (SI = 0.207, Figure 11a). However, the confidence intervals were wide (around 0.3 or higher) and overlapping.

**Figure 11.** Mean sensitivity indices and confidence intervals estimated via bootstrapping for NDVI in VBSA (**a**). Kolmogorov–Smirnov (KS) statistics for NDVI with regard to input factors in PAWN (**b**).

#### *3.5. Dense Mangroves—Normal Input Probability Distributions*

Compared with the DEN\_UNI case (Figure 5d,i), L2T(3) together with fCv(3) still had a dominant impact in NIR and SWIR for DEN\_NRM, but L2T(3) was less influential in some NIR bands in VBSA (Figure 5e,j). In addition, LIDFa(3) also had noticeable influences in NIR. WIDFa(3) became less influential compared with that for DEN\_UNI. Similar to SPS\_NRM, only a limited number of factors had positive SI values in VBSA for multiple VIs such as EVI, MNDWI and WFI (Figure 6e) and the confidence intervals were generally wide (0.25 or higher) in VBSA.

LIDFa(3) was investigated in this section as other influential factors such as fCv(3) and PAI(3) had been covered in previous sections. As shown in Figure 6e,j, NIR bands were the most sensitive to LIDFa(3). Therefore, the KS statistics and confidence intervals for S2A-B8 in PAWN were presented in Figure 12. fCv(3), L2T(3) and LIDFa(3) had high SI values and were distinct from other factors. The KS plots of fCv(3) and L2T(3) in Figure 12b were similar to those in Figures 8b and 9b, indicating that the output CDF diverged from the unconditional CDF when the fCv(3) and L2T(3) values increased or decreased. The KS plot of LIDFa(3) presented a similar pattern to that of fCv(3). Moreover, the KS values were higher when LIDFa(3) moved towards −0.6 (*θ*<sup>l</sup> = 67◦) than they were when LIDFa(3) was close to 0 (*θ*<sup>l</sup> = 45◦), which implied that canopy reflectance in NIR changed more significantly if leaves were vertically distributed.

**Figure 12.** Mean sensitivity indices and confidence intervals estimated via bootstrapping for S2A-B8 (**a**) and Kolmogorov–Smirnov (KS) statistics for S2A-B8 with regard to input factors (**b**) in PAWN.

#### *3.6. General Scenario with Correlated PAI(3) and fCv(3)*

If PAI(3) and fCv(3) were correlated for a general scenario, both were recognised as influential factors in PAWN. Besides, similar NSIs of canopy reflectance could be obtained with the GEN\_UNI with random PAI(3) and fCv(3) values (Figures 5f and 13b). In contrast, VBSA only identified PAI(3) as being influential in NIR and SWIR while fCv(3) had virtually no impact (Figures 5a and 13a). The computed NSIs in VBSA in Figure 13a might be meaningless because there were correlated inputs. Thus, no further interpretation was given. This plot demonstrated that the correlated biophysical properties of mangroves changed the calculated sensitivity indices in VBSA.

**Figure 13.** Stacked bar plots for normalised sensitivity indices of reflectance spectra in VBSA (**a**) and PAWN (**b**) for a general canopy with correlated PAI(3) and fCv(3). Each *x* value corresponds to an output, e.g., reflectance or VI, and each colour represents an input factor. The SI value is represented by the height of a bar and marked in the *y*-axis.

#### *3.7. A Brief Summary*

The results demonstrated the influences of SA methods and input configurations on the calculated SIs for reflectance of mangroves and the importance of examining the confidence intervals of calculated SIs, which might be ignored in previous SA studies on the RS of vegetation. The calculated NSI and SI values highlighted that canopy reflectance of mangroves was sensitive to fCv(3) and L2T(3) and PAI(3). The influences of *H*<sup>w</sup> for a sparse mangrove canopy and inclination distributions of plant material and *C*ab for a dense canopy might also be noteworthy. VIs with SWIR bands such as MNDWI and RGVI also had potential for mapping the fCv(3) and PAI(3) of mangroves and traditional VIs like EVI. Several VIs such as LSWI, MDI and WFI were sensitive to the L2T(3) of mangroves and might be helpful for estimating the LAI of mangroves from multispectral satellite images.

#### **4. Discussion**

#### *4.1. Global Sensitivity Analysis Methods and Interpretations of the Results*

The NSIs for canopy reflectance spectra of vegetation, as shown in Figure 5a–e, were frequently adopted in VBSA [19,20,22,47]. In addition, NSIs for VIs were also applied to evaluate the sensitivities of selected VIs to specific biophysical and biochemical properties of vegetation [22,47]. Stacked NSI plots for reflectance spectra were intuitive and straightforward for identifying the most influential factors and corresponding wavelength regions. However, SI values of various VIs could have distinct scales and only a limited number of factors might have non-negative SI values in VBSA, such as EVI and SAVI in Figure 6c. Therefore, a factor could obtain a high and dominant NSI value for a VI, e.g., the NSI of fCv(3) was 0.641 for EVI, while its absolute SI value was only 0.181 and not the highest compared with other VIs. Hence, the NSI values might mislead us to conclude that EVI was the most sensitive to fCv(3) in this case. Moreover, if a VI was only sensitive to a limited number of factors, the confidence intervals of these factors might be large and overlapping, and, thus, the results were not robust.

With regards to the output PDFs in VBSA, the results only showed that input distributions could affect output distributions but did not reveal any promising links between output PDFs and the confidence intervals of calculated SIs. A highly skewed output PDF such as RGVI in Figure 7y for SPS\_UNI did not necessarily mean the SIs in VBSA had wide

or overlapping confidence intervals. On the other hand, a normal-like output PDF such as S2A-B3 in Figure 7c for SPS\_UNI did not guarantee robust SIs in VBSA (Figure 9a). In addition, different numbers of input samples in VBSA did not significantly change the output PDFs as long as they were large enough but could result in divergent SI values in this study, especially for normally distributed input samples. Comparatively, the confidence intervals of calculated SIs for S2A band reflectance and VIs were usually narrow in PAWN with 2000 samples for the unconditional CDF and 200 samples per conditional point by 20 conditional points for conditional CDFs. The KS plots could be used to identify the value ranges of inputs where conditional output CDFs diverged from the unconditional CDF, which implied that there were noticeable changes in the output values.

#### *4.2. Differences between Sparse and Dense Mangrove Canopies*

The results also demonstrated that input ranges and distributions could affect the computed SIs, particularly in VBSA. For the GEN\_UNI case that PAI varied from 0 to 6, fCv(3), L2T(3) and PAI(3) were recognised as the most influential factors (Figure 5a,f). This partly coincided with the results of [20] that fCv and PAI were among the most influential factors. For the SPS\_UNI case, the impact of fCv(3) became slightly less dominant in VBSA while the impact of PAI(3) increased (Figure 5a,b). The changes in the impact of PAI(3) and fCv(3) in PAWN were not significant (Figure 5f,g). Additionally, increased influences of *H*w and fCv(1) in visible and NIR bands were identified by both methods. If the inputs were turned into normal distributions (SPS\_NRM), the computed NSIs of reflectance spectra also changed remarkably. As shown in Figure 5c,h, PAI(3) became less influential while the impact of fCv(3) and L2T(3) increased in NIR and SWIR, respectively.

If the canopy was dense with uniformly distributed inputs (DEN\_UNI), the NSI values of fCv(3) and PAI(3) decreased, particularly for PAI(3) (Figure 5d,i). Factors regarding the understory and background such as PAI(1) and *H*<sup>w</sup> had negligible influence now. In the meanwhile, biophysical properties related to woody material, such as L2T(3) and WIDFa(3), had a dominant impact in SWIR. The impact of LIDFa(3) also grew remarkably in NIR. In addition, the influences of leaf parameters, including *C*ab, *C*<sup>w</sup> and *C*m, increased by various degrees in different spectral regions. For the DEN\_NRM case, the NSI patterns in PAWN (Figure 5j) were similar but slightly different from Figure 5i for DEN\_UNI. The impact of fCv(3), L2T(3) and LIDFa(3) increased while that of WIDFa(3) decreased. In VBSA, fCv(3) and L2T(3) became the most influential factors, and the influences of LIDFa(3) and WIDFa(3) both decreased (Figure 5e).

In summary, PAI(3) had a larger impact under the sparse canopy while the influences of leaf parameters and inclination angles increased under a dense canopy. This agreed with the results in [19,22]. The new factors fCv(3) and L2T(3) were among the most influential factors under all examined scenarios except for the correlated case where fCv(3) was not identified as influential in VBSA. *H*<sup>w</sup> mainly affected the reflectance of sparse canopies, but it was worth noting that the water body was not allowed to immerse the crown in this study. Otherwise, the infrared reflectance of submerged mangroves would reduce dramatically, which had been used to identify them from satellite images [4,48].

#### *4.3. Potential Limitations and Suggestions*

The structure of real mangrove forests was more complicated than the model could simulate. Simplifications in the model might result in overestimated contributions from fCv(3) because it was applied to weight the reflected radiation from mangroves directly. Additionally, the woody material was assumed to be randomly distributed and to share the same fractional cover with leaves in the adopted model [28]. Its reflectance spectra were averaged from field measurements over the lower parts (lower than 2 m) of stems [28], which might be different from the reflectance of branches and shoots in the crowns. Consequently, the contributions of L2T(3) and its SIs could also be overestimated.

Although there were potential divergences in the computed SI values of L2T(3), the impact of woody material on canopy reflectance were demonstrated [5,49]. In addition, caution for the inclination distributions of leaves and wood was still needed because they could be an obstacle when mapping the LAI of mangroves from satellite images. The leaf area ratio could be measured by classifying the point clouds acquired by laser scanning [50,51]. Average inclination angles could also be calculated from point clouds [52,53] or digital photos [54]. However, convenient and operational methods were still lacking [55], particularly when they needed to be measured together with PAI and fCv at many of the plots to be used as calibration or validation data.

For future studies, it might be a good start to map the fCv of mangroves together with PAI because field data of fCv were relatively easy to collect, e.g., by taking hemispherical photos or upward photos [56]. Moreover, the data could be used to assess the potential correlation between PAI and fCv, and further the impact of fCv and background water on PAI mapping. According to the GSA results in this study, multiple VIs may be suitable for mapping PAI and fCv of mangroves from S2A images. Besides traditional VIs such as EVI and NDVI, other VIs such as MNDWI and RGVI were also sensitive to fCv and PAI. The new VIs for mangroves or aquatic plants such as MDI, WFI and RGVI had a common characteristic that SWIR bands were incorporated. This study showed that the SWIR bands and these VIs could be helpful for mapping mangroves, but the optimal VI for mapping a biophysical property from a specific dataset might vary.

#### **5. Conclusions**

Based on a canopy reflectance model of mangroves, this study demonstrated that GSA methods (variance-based or density-based), input ranges and probability distributions (uniform or normal) could affect the computed sensitivity indices under the examined mangrove scenarios. Briefly, fCv(3) and L2T(3) were among the influential factors for the infrared reflectance of mangroves under the examined scenarios. PAI(3) was also influential for a sparse canopy but became less influential for a dense canopy. In contrast, inclination distributions of plant material and leaf parameters, e.g., leaf water content, could become more influential in infrared bands for a dense canopy. Moreover, the influence of water depth was noteworthy for a sparse canopy and maybe other scenarios if the water body could immerse the crown.

Since the results and conclusions can be different if a specific model, GSA method and input configuration are adopted, it may be essential to perform a tailored GSA according to the study area and available data. Based on the results, it is recommended that attention should be paid to the L2T and fCv as they may affect estimating the LAI or PAI of mangroves. A complete field dataset that includes the factor to be mapped such as LAI and if possible, other influential factors such as fCv, L2T ratio and inclination angles would be beneficial to further analyses and evaluations of their impact. This, in turn, could contribute to the development of a protocol on field data collection for mapping mangrove from remote sensing images. Considering the challenges in field data collection, it may be a good start to collect and map PAI and fCv of mangroves. For VI based mapping methods, more choices such as MNDWI and RGVI deserve an attempt if SWIR bands are available, but tests are still needed to determine the optimal VI for mapping a specific property based on the available field data and RS images.

**Author Contributions:** Conceptualization, C.N. and S.P.; Data curation, C.N.; Formal analysis, C.N.; Funding acquisition, C.N. and S.P.; Methodology, C.N.; Project administration, C.N. and S.P.; Resources, C.N. and S.P.; Software, C.N.; Supervision, S.P. and C.R.; Writing—original draft, C.N.; Writing—review and editing, S.P. and C.R. All authors have read and agreed to the published version of the manuscript.

**Funding:** Chunyue Niu was supported by an Australian Government Research Training Program (RTP) Scholarship and Remote Sensing Research Centre, School of Earth and Environmental Sciences, The University of Queensland.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** We thank Francesca Pianosi, Fanny Sarrazin and Thorsten Wagener at the University of Bristol for sharing the SAFE toolbox. We are grateful to the anonymous reviewers for their valuable comments.

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

#### **References**


## *Article* **Mangrove Forest Cover and Phenology with Landsat Dense Time Series in Central Queensland, Australia**

**Debbie A. Chamberlain 1,2,\*, Stuart R. Phinn <sup>2</sup> and Hugh P. Possingham 1,3**


**Abstract:** Wetlands are one of the most biologically productive ecosystems. Wetland ecosystem services, ranging from provision of food security to climate change mitigation, are enormous, far outweighing those of dryland ecosystems per hectare. However, land use change and water regulation infrastructure have reduced connectivity in many river systems and with floodplain and estuarine wetlands. Mangrove forests are critical communities for carbon uptake and storage, pollution control and detoxification, and regulation of natural hazards. Although the clearing of mangroves in Australia is strictly regulated, Great Barrier Reef catchments have suffered landscape modifications and hydrological alterations that can kill mangroves. We used remote sensing datasets to investigate land cover change and both intra- and inter-annual seasonality in mangrove forests in a large estuarine region of Central Queensland, Australia, which encompasses a national park and Ramsar Wetland, and is adjacent to the Great Barrier Reef World Heritage site. We built a time series using spectral, auxiliary, and phenology variables with Landsat surface reflectance products, accessed in Google Earth Engine. Two land cover classes were generated (mangrove versus non-mangrove) in a Random Forest classification. Mangroves decreased by 1480 hectares (−2.31%) from 2009 to 2019. The overall classification accuracies and Kappa coefficient for 2008–2010 and 2018–2020 land cover maps were 95% and 95%, respectively. Using an NDVI-based time series we examined intra- and inter-annual seasonality with linear and harmonic regression models, and second with TIMESAT metrics of mangrove forests in three sections of our study region. Our findings suggest a relationship between mangrove growth phenology along with precipitation anomalies and severe tropical cyclone occurrence over the time series. The detection of responses to extreme events is important to improve understanding of the connections between climate, extreme weather events, and biodiversity in estuarine and mangrove ecosystems.

**Keywords:** Landsat; mangrove forests; time series; Google Earth Engine; random forests; phenology; TIMESAT; climate; monitoring; Great Barrier Reef

#### **1. Introduction**

Coastal and estuarine ecosystems are valued for their ecosystem services, particularly coastal protection, as a nursery for fish, and carbon sequestration [1,2]. Mangrove forests are an important part of tropical estuarine ecosystems and considered integral to the emerging blue carbon economy. Blue carbon consists of carbon that is stored, sequestered, or released from coastal vegetation ecosystems [3]. Despite this, human interference leading to widespread degradation and deforestation has caused a decline in mangrove cover and biomass. Carbon emissions from global mangrove loss is estimated to be 2391 Tg CO2 eq by the end of the century [4]. For reference, emissions from the global transportation sector are projected to rise to 11,900 Tg CO2 eq by 2100 [5]. Although the rate of mangrove loss has decreased substantially since the 1990s, from ~2% to <0.4% per

**Citation:** Chamberlain, D.A.; Phinn, S.R.; Possingham, H.P. Mangrove Forest Cover and Phenology with Landsat Dense Time Series in Central Queensland, Australia. *Remote Sens.* **2021**, *13*, 3032. https://doi.org/ 10.3390/rs13153032

Academic Editor: Brigitte Leblon

Received: 8 May 2021 Accepted: 28 July 2021 Published: 2 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

year [6] increasing human pressure along the coast has synergistic effects with climate change and threatens the distribution and abundance of mangroves. As coastal forests with unique adaptations to saline conditions, mangroves form a characteristic vegetation zone along sheltered bays, tidal inlets and estuaries in the tropics and subtropics [7]. The major threats to mangroves are climate change, clearing (through urban areas, ports, aquaculture, and industry development), timber harvesting, dieback, changes in hydrology (e.g., the restriction or alteration of flows), and pollution [8–10]. Understanding mangrove ecosystems and mapping their extent is critical to meeting the UN sustainable development indicator 6.6.1: "change in the extent of water-related ecosystems over time". Indicator 6.6.1 is used in determining progress toward meeting Sustainable Development Goal 6 (SDG 6), which is to "ensure availability and sustainable management of water and sanitation for all" [11]. Forest degradation mapping at the canopy scale supports SDG 6, and is primarily based on time series of vegetation indices obtained from optical sensors on satellites, e.g., medium spatial resolution Landsat satellites [12].

A key requirement for understanding the impact of human activities on climate, ecosystem function, and biodiversity is assessing the rate of change in tropical forest cover [13,14]. Many studies have analysed change dynamics in forest and wetland ecosystems through remote sensing data over time [15,16]; however, image sequences are often incomplete and of limited temporal range [17,18]. In tropical regions, many observations each year are required to obtain a complete cloud-free set of images for a region of interest [13,19]. Landsat-based classification with supervised machine learning techniques for land cover, along with spectral transformations using Vegetation Indices (VIs), are an efficient method for accurate mapping and monitoring of tropical forests and mangroves [20–22]. Random Forests (RF) are an ensemble learning technique for land cover classification, robust to outliers and noise, and computationally undemanding than other classification methods (e.g., boosting, support vector machines) for practical applications with a lot of data [23–25]. Given the reported high accuracy of RF in land cover classification of forests and wetlands, RF have been employed to map the extent of mangrove forests in both the sub-tropics and tropics [26–28].

The most commonly used VI is the normalized difference vegetation index (NDVI) [29], a broadband green-sensitive (photosynthetically active) index for detecting seasonal and inter-annual variations in canopy vigour and growth in relation to climate parameters and relating these variations to the capacity of the canopy to photosynthesize [30]. In accordance, phenological variations are changes in the rate of growth as indicated by changes in photosynthesis that can be detected by satellites through spectral transformations of VIs. Indeed, the variation from season to season in the satellite spectral data as measured by remote sensing and vegetation phenology, is an essential approach for understanding the role of vegetation in climate change [31]. Xiao et al. [32] challenged the idea that habitats in tropical environments do not show variation from season to season and they identified changes in growth rate driven by variation in rainfall. With dense time series, the suitability of phenology-based smoothing techniques such as harmonic regression and Savitsky-Golay have recently been demonstrated to map tropical forest disturbance and degradation [33,34]. Furthermore, to capture the temporal variation of vegetation growth cycles the localised climate at each time series location should be extracted. Localised precipitation patterns will influence the phenophases, and if the spatial distribution of weather stations is insufficient, these subtle variations in climate will not be recorded [35]. Finally, it has been shown that under a comparable climatic environment, genotypes of species distributed across a study region may vary enough to express a different phenophase [36].

Few studies have used VIs as phenological variables in classification approaches to maximise spectral contrast among vegetative communities and discriminate wetlands and mangroves from other vegetation types (e.g., grasses, crops) [16,37]. Li et al. [38] investigated wetland dynamics with greenness trends under intense land cover change on the coast of China. Lamb et al. [39] mapped estuarine emergent wetlands based on seasonal structural change, and [40] identified salt marsh seasonal trend decomposition in the Chesapeake and Delaware Bays, USA. Wu et al. [16] extracted phenological differences in the wetland vegetation of Chongming Island, China with time series harmonic regression. As these analyses were computationally demanding, the researchers used time series of the Landsat archive linked with wetland vegetation phenology features as a Google Earth Engine (GEE)-managed pipeline. The GEE platform enables processing and classification applications for large-scale geospatial mapping and monitoring of land cover and land use change [22]. Nonetheless, change dynamics and degradation in mangrove forests integrating phenology-based dense time series has rarely been investigated with remote sensing approaches [15]. In this paper, we use the term 'phenology' in a general context to describe the cyclic vegetation activity of mangroves over time. Here, we apply machine learning classification procedures to a time series of mangrove cover and examine climate drivers of mangrove growth dynamics across Central Queensland, Australia. To our knowledge, this is the first study to examine mangrove phenology across different sites on the east coast of Australia. Harmonic regression captures cyclic behaviour on an intraand inter-annual basis, and the Savitzky-Golay filter captures subtle dynamics during the vegetative growing season. We used both approaches to fully examine the seasonal variability in mangrove forests in our study region. The study region is a large estuarine area adjacent to the Great Barrier Reef, Australia. Although management of Great Barrier Reef catchments is important to maintain the ecological integrity of the reef, monitoring of estuarine areas receives little attention. Moreover, monitoring of Ramsar Wetland sites on the east coast is absent. The objectives of the study are: (1) to quantify how coastal mangrove forests have spatially and temporally changed in a decadal period (2009–2019) within estuarine catchments of the Great Barrier Reef; (2) to investigate the seasonality in tropical mangrove forests through a dense time series of satellite images, exploring the inter-annual and seasonal responses and variations of mangrove phenology to 21 years of rainfall variability; (3) to inform regional coastal planning, conservation efforts, and policy-makers.

#### **2. Materials and Methods**

#### *2.1. Study Area*

Our study area is located within the northeast coast drainage division of the Central Queensland coast encompassing the Mackay Whitsunday Natural Resource Region and part of the Isaac Region, Central Queensland, Australia (Figure 1). Cape Palmerston National Park is positioned in the Ince Bay Receiving Waters adjacent to the World Heritage listed Great Barrier Reef. The Shoalwater and Corio Bays Area Ramsar site is located south of Cape Palmerston National Park. A large proportion of the marine waters in the Ramsar site are included in marine parks (Commonwealth and Queensland), including the Great Barrier Reef Marine Park (Commonwealth) and Great Barrier Reef Coast Marine Park (State) [41]. The Ramsar wetland supports a broad range of natural values including nationally/internationally threatened wetland species, significant species diversity and large populations of waterbirds, green turtles, dugong, and fish that use the site for vital life history functions such as roosting, nesting, feeding, and breeding. Cape Palmerston National Park is listed as a category II protected area on the International Union for Conservation of Nature (IUCN) World Database on Protected Areas and covers 7200 hectares [42]. The Shoalwater and Corio Bays Area Ramsar site covers 239,100 hectares [43].

The primary intensive land use in the region is the cultivation of sugar cane, making up 18% of the catchment area, with Mackay being the largest sugar-producing region in Australia [44]. Grazing is also an important land use, accounting for 42% [45]. The region's estuaries directly support several commercial fisheries, e.g., East Coast Inshore Fin Fish Fishery, East Coast Otter Trawl Fishery, and Line Fishery (Reef) [46]. Additionally, recreational fishing is a considerable activity in the region, with 25% of the population participating in fishing for recreation, far greater than the state average of 15% [47]. Mangroves and associated communities cover 87,994 hectares of tidal land in the region, with nine wetland areas recognized as nationally important [48], and one Ramsar wetland

recognized as internationally important. The Ramsar wetland site is subtropical with species from both temperate and tropical regions. Within the high rainfall areas of the Central Queensland coast bioregion, estuarine wetlands are about equally dominated by saltpan and samphire flats along the high intertidal area; yellow and orange mangroves (*Ceriops tagal* and *Bruguiera* spp.) dominate along the mid-intertidal area; and the stilted mangrove (*Rhizophora stylosa*) dominates in the lower intertidal area [48]. The river mangrove *Aegicera corniculatum* occupies an intermediate-upstream estuarine position while the ubiquitous grey mangrove *Avicennia marina* inhabits a wide range of tidal inundation levels [49]. Twenty-three tree and shrub species of mangroves are present [50]. Mean maximum temperature is between 22.7 ◦C in July and 30.4 ◦C in January, while the mean minimum temperature is between 11.4 ◦C in July and 23.2 ◦C in January. The dry season is May-November, the wet season occurs December-April. Mean annual rainfall is 1528.4 mm, while the wettest month (344.8 mm) is recorded in February. The study area is located between latitude 20◦22 –22◦58 S and longitude 147◦19 –150◦46 E (Figure 1).

**Figure 1.** Study site with weather stations at Lethebrook, Koumala, and St. Lawrence-Bowen to Shoalwater and Corio Bays Area Ramsar Wetland, Central Queensland-Landsat fractional cover image visualised by using the Band 1-red (bare ground and rock), Band 2-green (vegetation), Band 3-blue (non-green vegetation indicative of drier habitats with less vegetative cover), provided by [51].

#### *2.2. Image Classification*

Medium-resolution satellite imagery is suitable for mapping mangrove and wetland areas on a regional scale [52]. There are two reasons for selecting Landsat imagery: (1) it is acquired at regular intervals, and (2) it is freely available from USGS. Data processing and analysis was performed in Google Earth Engine (GEE). We used a Landsat time series to understand how mangroves have changed over time, identify areas of loss/gain and understand patterns of change in our study region. Landsat data were derived from the United States Geological Survey (USGS) Landsat 7 and 8 Collections, atmospherically corrected surface reflectance products [53]. In regions that are severely cloud-affected, such as the tropical coastal zone, it is necessary to accumulate enough Landsat data to create a cloud-free composite, therefore we chose three-year decadal windows, 2008–2010 and 2018–2020. We created a composite on a per pixel per band basis using the median reducer in GEE [27]. Topographic variables were included in the GEE image stack as a 1 arc-second (approximately 30 m) digital elevation model, the Shuttle Radar Topography Mission (SRTM), sourced from the National Aeronautics and Space Administration Agency (NASA) Jet Propulsion Lab [54]. The SRTM provides elevation and slope data subset to include only low elevation coastal areas where mangroves may be present. We masked to elevations in our study region that are less than 65 m. above mean sea level. Two distinct data sources served as reference datasets: Global Mangrove Watch (GMW) Project [55] providing the 2009 reference data; and Geoscience Australia's Digital Earth Australia platform [56] providing the 2019 reference data. For cloud/shadow removal we filtered the Landsat composite using the quality assessment (QA) band [57] and the GEE median reducer. We calculated four VIs [58,59]: the Simple Ratio (SR) indicates healthy vegetation; the Normalized Difference Vegetation Index (NDVI), the most widely used VI and useful in understanding vegetation density and assessing changes in plant health; the Modified Normalized Difference Moisture Index (MNDWI) assesses vegetation water content; Green Chlorophyll Vegetation Index (GCVI) estimates green leaf biomass. In GEE, we applied masks to areas of low elevation and high NDVI and MNDWI. Additional masking allows us to focus on areas that are more likely to have mangroves. The auxiliary variables in the analysis were the topographic variables derived from the digital elevation model. The phenology variable was the median NDVI calculated from the time series in the three-year decadal windows, 2008–2010 and 2018–2020. The phenology variable was included to separate mangrove trees from other vegetation types (e.g., grass and crops) based on differences in seasonality [60].

Since this study seeks to evaluate predictor variables that tell us the extent of mangrove forest, only two land cover classes were generated (mangrove versus non-mangrove that includes bare soil, sand beach, water, crops, pasture, and native vegetation). We sought a supervised classification approach that is robust to normal distribution departures, can overcome limitations such as overfitting, and uses different subsets of the same training dataset [61]. Machine learning allows us to use samples of areas with and without mangroves to detect mangrove forests across a region. RF is a supervised, nonparametric, and ensemble tree-based machine learning algorithm, which has increasing use in mangrove mapping with Landsat images [62]. RF has been found to outperform comparative classification models in projection accuracy and computation cost [63,64]. Based on the predictors (Landsat bands), the trees will vote for each pixel to detect mangrove vs. non-mangrove with the most supported value assigned to each pixel. We used 100 decision trees with five predictors including the vegetation indices NDVI, MNDWI, SR, and GCVI. In GEE, the RF classifier chooses five predictors based on a random selection from seven bands in Landsat 8: 'B5', 'B6', 'B4', and the VIs. In the same way, the RF classifier chooses five predictors based on a random selection from seven bands in Landsat 7: 'B4', 'B5', 'B3', and the VIs. The training samples were collected by creating geometries of the location of mangroves using the VIs and elevation masks with Landsat false-colour composites in GEE. The total number of sample points (pixels) was 19,500. In our analysis, 80% of the sample points were randomly split and used to train the classifier. The decision trees were fully grown without pruning using a sample (with replacements) of 20% of the training data, also randomly split, to test the accuracy and validate the RF classifier [65,66]. For accuracy assessment a confusion matrix was calculated based on the sample data and classification result giving the evaluation metrics, user's and producer's accuracy, overall accuracy, and Kappa coefficient. Moreover, the mangrove forest maps were visually inspected and compared with high-resolution mangrove canopy cover maps from Digital Earth Australia to evaluate the performance of our method [56]. Considering that the Digital Earth Australia maps quantify the change in the extent of mangrove forest over the period 1987 to 2016 at an annual, national scale, there is concurrence with the classification maps in our study that demonstrate a regional and local approach.

#### *2.3. Mangrove Phenology*

All available images (i.e., 2012 to 2020) for our study region (Bowen to Shoalwater and Corio Bays Area Ramsar Wetland) were acquired from the Landsat 8 surface reflectance data with a spatial resolution of 30 m and temporal resolution of 16 days and pre-processed in GEE using functions including masking clouds/cloud shadows. Images were stacked together to build a time series to be used in the models.

We used linear and harmonic regression models to analyse the seasonal parameter trends in our time series. First, we fitted an ordinary least squares model to the time series of all Landsat 8 statistical coefficients with the linear regression reducer in GEE. Second, we detrended the data to reveal underlying fluctuations in the time series. Third, we fitted Fourier-style harmonic regression equations including a trend term to each pixel and spectral component to identify phenological behaviours [37]. Harmonic regression is a mathematical technique largely focused on NDVI that is used to decompose a complex, static signal into a series of individual cosine waves (terms) and an additive term [31]. The time series is approximated as a trigonometric polynomial, with coefficients added to the harmonic model to get the amplitude (A) and phase (ϕ):

$$p\_t = \beta\_0 + \beta\_1 t + \text{A } \cos(2\pi \omega t - \varphi) + e\_t = \beta\_0 + \beta\_1 t + \beta\_2 \cos(2\pi \omega t) + \beta\_3 \sin(2\pi \omega t) + e\_t \tag{1}$$

where *pt* is the value of the harmonised NDVI within a pixel obtained from the satellite data at time *t*, β<sup>0</sup> is a constant component of the regression that designates the start of greenness at each pixel, β<sup>1</sup> is a linear coefficient (slope/linear trend), β<sup>2</sup> is the cosine term coefficient, β<sup>3</sup> is the sine term coefficient, *et* is a random error, and ω is the angular frequency. We set ω = 1 (one cycle per unit time for an annual cycle) to fit the model to the time series.

#### TIMESAT Analysis

The program TIMESAT is an approach for extracting phenological parameters from remote sensing data and was partly developed from the premise that changes in vegetation phenology may be an indication of climate change [67]. The method is based on nonlinear least squares fits to the upper envelope of the NDVI data [68]. First, we applied a median filter and interpolated the missing values. Second, we used the Savitzky-Golay model and seasonal trend decomposition for a yearly growing season in tropical mangrove forests. The Savitzky-Golay filter follows within-season variations and therefore captures subtle dynamics during the growing season [69]. The result is a smoothed curve adapted to the upper envelope of the NDVI values [68]. To analyse the dense time series (1999–2020), we first created a harmonise function between L5, L7, and L8 sensors. There are small, but potentially significant differences between the spectral characteristics of Landsat ETM+ and OLI, depending on the application. It is advisable to harmonise for a long time series that spans Landsat TM, ETM+, and OLI sensors [31]. Following the harmonise function, we added the NDVI band to the collection. We extracted 13 phenological metrics for each season (21 years) from the TIMESAT program with Landsat dense time series at our study region. TIMESAT provides phenology parameters such as: start of season (SOS) depicted as initiation of a consistent upward trend in the NDVI time series, or the beginning of measurable photosynthesis in the mangrove canopy; end of season (EOS) identified as the end of consistent downward trend in the time series, or the end of measurable photosynthesis in the mangrove canopy; season length (LOS) or duration of photosynthetic activity (time of SOS to EOS); peak NDVI value (the maximum recorded NDVI value); base value (the average of the minimum values left and right of the season); amplitude (the range between the maximum NDVI and the base level or the maximum increase in

canopy photosynthetic activity above the baseline); peak time (POS: time of the middle of the season, or time of maximum photosynthesis in the canopy); left derivative (rate of increase at the beginning of season or growth rate); right derivative (rate of decrease at the end of season or rate of senescence [70–72]. Mangrove productivity over the region was estimated using the amplitude, peak NDVI, left and right derivates, large integral, and small integral. The large integral is the function describing the season from the start to the end (an estimate of the total vegetation productivity in the annual cycle). The small integral is the difference between the function describing the season and the base level from season start to end. It is a measure of the productivity within the growing season [70,73]. Further, we used the Seasonal-Trend Decomposition by LOESS (STL) method to remove spikes and outliers in the time series and produce a smoother trendline [74].

Relationships between TIMESAT seasonality parameters were examined with linear regression models by the increasing/decreasing patterns (regression slope) and goodnessof-fit (adjusted coefficient of determination, Adjusted R2). To characterise the geographic distribution of environmental drivers of mangrove phenology, we applied climate data for 1999–2020 to the time series. Monthly rain gauge precipitation data interpolated from ground station measurements were obtained from the Australian Bureau of Meteorology online climate database (www.bom.gov.au, accessed on 1 March 2021). For the northern section of our study region, rainfall data were primarily gathered from adjacent Lethebrook Station (station number: 33041) with missing data supplemented from Bowen Station (station number: 33264). For the mid-section of our study region rainfall data were primarily gathered from the adjacent Koumala Station (station number: 33038) with missing data supplemented from Plane Creek Station (station number: 33059). For the southern section of our study region rainfall data were collected from two adjacent stations at St. Lawrence (station numbers: 33065 and 33210). As rainfall amount differs across the region, trendlines of the precipitation data were superimposed on the interpolated and filtered NDVI from the three locations across our study site: northern at Repulse Bay, mid-section at Rocky Dam Creek/Cape Palmerston National Park, and southern at Herbert Creek bordering Shoalwater and Corio Bays Area Ramsar Wetland. The peak of the season (POS) in relation to precipitation trends were compared between sites.

RF classification, linear and harmonic regression modelling and dense time series data processing was done with Google Earth Engine. The TIMESAT program was used for phenometric extraction (Figure 2). Statistical analysis was done with the R statistical environment [75] (r-project.org/).

**Figure 2.** Workflow used for land cover classification and phenological analysis.

#### **3. Results**

The aim of this research is to quantify how coastal mangrove forests have spatially and temporally changed in a decadal period (2009–2019) within estuarine catchments of the Great Barrier Reef and to explore the seasonality, climate responses, and variability of mangroves through a dense time series of satellite images in Central Queensland, Australia. Both biotic and abiotic factors constrain the phenological trajectory of a mangrove forest, e.g., regional climate, seawater and soil salinity, latitude, and local vegetation interactions, and these can vary across different locations and years. Sea surface temperatures in the Australian region have warmed by around 1 ◦C since 1910, with the Great Barrier Reef warming by 0.8 ◦C in the same period. Moreover, observations show that there has been an increase in the intensity of heavy rainfall events in Australia. (Figures S1 and S2 in the Supplementary Materials). Against the background of climate change, our study examines the phenodynamics of mangrove forests across different sites on the east coast of Australia.

#### *3.1. Image Classification*

The mangrove forest extent derived from the RF classification shows a decrease of 1480 hectares (−2.31%) over the decadal window, 2009–2019 (Table 1). In ArcGIS, we generated a Landsat 8 Operational Land Imager (OLI) true-colour mosaic of our study region, overlaying the mosaic with the mangrove change map from GEE (Figure 3). The overall classification accuracy and Kappa coefficient for the 2009 and 2019 land cover maps were 95% and 95%, respectively (Table 2). The high accuracies suggest that the RF classifier is an efficient model for projecting mangrove forest cover dynamics on both training and test datasets. The mangrove forest maps satisfactorily aligned with the high-resolution mangrove canopy cover continental-scale maps from Digital Earth Australia [56].

**Table 1.** Mangrove forest extent derived from Random Forest classification with Landsat NDVI images 2009 and 2019 (Bowen to Shoalwater and Corio Bays Area Ramsar Wetland, Central Queensland).


**Table 2.** Accuracy assessment from Random Forest classification 2009 and 2019 Landsat decadal time series analysis (Bowen to Shoalwater and Corio Bays Area Ramsar Wetland, Central Queensland).


**Figure 3.** Landsat Operational Land Imager (OLI) true-colour mosaic (images captured 16 July 2020 and 9 August 2020) overlaid with the change map of mangrove forest extent (decline) in red from the Random Forest classification and showing the three sites input to GEE for time series construction and subsequent phenology analysis: northern site, Repulse Bay; mid-section, Rocky Dam Creek/Cape Palmerston National Park; and southern site, Herbert Creek, study region is Bowen to Shoalwater and Corio Bays Area Ramsar Wetland.

#### *3.2. Mangrove Phenology*

The detrended time series (2012–2020) shows peaks and troughs in the time series but has missing data caused by cloud masking of pixels in GEE (Figure S3 in the Supplementary Materials). Nevertheless, there is a decrease in mean NDVI across the time series with strong seasonality. Values of the peak growth and photosynthetic activity range from 0.72 to 0.9 and values of the minimum photosynthetic activity range from 0.62 to 0.82. The predicted values of the mangrove harmonic regression model were compared with the actual mean NDVI values and show a general trend in the time series (Figure 4). The seasonal trend indicates that high values were attained towards the end of the growing season (generally in winter) and low values were realized in late spring (October to November with the exception of 2018 where the lowest value was in summer, i.e., December, and 2019 where the lowest value was in early spring, i.e., September, Figure 5). We suggest that under tropical conditions, frequent cloud cover and high temperatures in spring and summer may limit the ability of mangroves to photosynthesize hence greater NDVI values occur in winter. This outcome agrees with field studies conducted by [49] with *Avicennia marina* in tropical Australia. We identified that tropical mangroves on the east coast of Australia display a disproportional bimodal seasonality with two periods of high mean NDVI and two periods of low NDVI values annually through the Landsat time series (1999–2020) (Figures 6 and 7).

**Figure 4.** Harmonic regression model original and fitted values of all available (2012–2020) Landsat 8 imagery with mangrove NDVI. The red line is the fitted NDVI values and the blue dotted line is the original NDVI values for mangrove forest in the mid-section of our study region. The site is located directly adjacent to Koumala Station, Figure 1, Bowen to Shoalwater and Corio Bays Area Ramsar Wetland, Central Queensland.

**Figure 5.** Mean NDVI temporal pattern across the dense Landsat time series (1999–2020) for mangrove forest in the mid-section of our study region. The site is located directly adjacent to Koumala Station, Figure 1, Bowen to Shoalwater and Corio Bays Area Ramsar Wetland, Central Queensland.

**Figure 6.** Mangrove forest intra-annual smoothing phenology with the Savitzky-Golay filter in TIMESAT using Landsat NDVI dense time-series data 1999–2020, blue line indicates the decadal NDVI data, brown line indicates the smoothed time series fitted data, brown dots display the start and end of the growing season for mangrove forest in the mid-section of our study region, Figure 1, Bowen to Shoalwater and Corio Bays Area Ramsar Wetland, Central Queensland.

**Figure 7.** Mangrove forest seasonal trend decomposition in TIMESAT with Landsat NDVI data and dense time-series 1999–2020, black line indicates the seasonal trend, grey dots display the start and end of the growing season for mangrove forest in the mid-section our study region, Figure 1, Bowen to Shoalwater and Corio Bays Area Ramsar Wetland, Central Queensland.

#### TIMESAT Analysis

In the TIMESAT program we used the Savitzky-Golay algorithm and extracted 13 attributes describing the spectro-temporal profile of mangroves for each season (Table 3). The filtering process preserved the trend of the original NDVI data and the fitted curve depicted the seasonal variation in forest characteristics; it, therefore, provides reliable data for the extraction of phenological information. For the mid-section of our study region (Rocky Dam Creek/Cape Palmerston National Park, Figure 3) we found the SOS to be December, EOS to be September/October, and the LOS to be 10 months. POS varied through the timeseries (February to July): 1999–2003 in June/July; 2004 in February; 2005–2006 in May; 2007–2008 in June/July; 2009–2010 in May/June; 2011–2012 in May; 2013 in April; 2014–2020 in June/July (Figure 6). There is a marked upward movement in the seasonal trend decomposition of NDVI data from the summer of 2009 to 2020 (Figure 6). Notwithstanding the length of the season gradually becoming shorter, the long integral (total productivity in the annual cycle) increased (Figure S4 in the Supplementary Materials). The trend in the peak NDVI value (maximum recorded NDVI value in the time series) increased from the summer of 2009 (Figure 7 and Figure S5 in the Supplementary Materials). In evergreen areas the small integral (productivity within the growing season) may be small even if the total vegetation production is large. The small integral and amplitude (difference between the maximal value and the base level) both declined from 2009 in our study region due to a positive change in the base value increasing the amount of green vegetation for this period (Figures S6 and S7 in the Supplementary Materials). We found a significantly high correlation (*p* < 0.001) between the phenology metrics from TIMESAT through the time series particularly season length and duration of the peak value (Pearson correlation coefficient r value = −0.91, Adjusted R<sup>2</sup> = 0.8182), and season (year) and season length (Pearson correlation coefficient r value = −0.91, Adjusted R2 = 0.8179) indicating a dynamic trend in the cyclical behaviour of tropical mangrove forests (Table 4).

Except for the northern region, Central Queensland is largely decoupled from southern monsoonal influences and exhibits highly variable phenology that is rainfall pulse driven. Monthly rainfall was 0–800 mm at our study site, Rocky Dam Creek/Cape Palmerston National Park (Figure 3). Australia's millennium drought of 2001–2007 is reflected in the rainfall pattern and the reduced mean NDVI trendline (Figure 8). Considerable decreases in NDVI occurred in 2001, 2002, 2004, and 2008 that are consistent with the trend of monthly rainfall. In addition to the general trend, phenological characteristics and rainfall patterns correspond with major climate events such as the EL Niño-Southern Oscillation (ENSO). The increase in precipitation from the summer of 2009 in the mid-section of our study region produced a concomitant increase in photosynthetically active growth and biomass production in mangrove forests at this site. The extreme and extended rainfall event from La Niño that occurred over Queensland in the summer of 2010/2011 was reflected in the net seasonal photosynthetic rates for this period. Further, a strong La Niño event is apparent in the Southern Oscillation Index (SOI, the difference in surface air pressure measured between Tahiti and Darwin, northern Australia) that spikes in the 2010/11 summer (Figure S8 in Supplementary Materials). In addition, a severe tropical cyclone (Tropical Cyclone Yasi, Category 5) brought an intense rainfall event in February 2011 to the entire region [76]. A relatively high rainfall pattern continued until mid-Autumn of 2017 when a severe tropical cyclone (Tropical Cyclone Debbie, Category 4) impacted the coast and destroyed coastal vegetation [77]. Recovery from the cyclone is evident in 2018 when the mean NDVI reached 0.9 at Rocky Dam Creek/Cape Palmerston National Park (Figure 8).


 of our study region, Startt.:

photosynthesis

 in the

**Table 3.**

start of season is the beginning of measurable

Savitsky-Golay

 filter with

phenometrics

 extracted from TIMESAT for the Landsat dense time series of mangrove forest in the mid-section

photosynthesis

 in the mangrove canopy (month of the year); Endt.: end of season is the end of measurable


**Table 4.** Response and explanatory variables used in the linear regression models with Pearson correlation coefficient r values and adjusted R2 at significantly high values (*p* < 0.001).

**Figure 8.** Interpolated, filtered NDVI and rainfall trendlines for the mid-section of our study region, Rocky Dam Creek/Cape Palmerston National Park, rainfall data gathered from Koumala Station (station number: 33038).

> In comparison to the mid-section, the northern section of our study region, Repulse Bay, displays a generally higher level of mean NDVI and less temporal variability (Figure 3 and Figure S9 in the Supplementary Materials). Monthly rainfall was 0–900 mm across the time series. The POS occurred in diverse months (May to August), exemplifying a different

climate pattern at this site. The northern site exhibited a NDVI decrease in 2002 similar to our mid-section site; however, the decrease occurred in the winter of 2002 (at the northern site it was early summer) and was less pronounced (Figure 9). From 2003–2020 a steady increase in NDVI activity was apparent. High rainfall peaks occurred from 2007–2012, followed by 6 years of moderate rainfall, escalating to a monsoonal trough event in January 2019, then falling to moderate rainfall conditions in 2020. A tropical cyclone (Tropical Cyclone Jasmine, Category 2) affected the northern section in February 2012. A substantial decrease in NDVI was apparent in the winter of 2017 that was not reflected in the rainfall regime but possibly occurred from the influence of cloud cover.

**Figure 9.** Interpolated, filtered NDVI and rainfall trendlines for the northern section of our study region, Repulse Bay, rainfall data gathered from Lethebrook Station (station number: 33041).

> The southern section of our study region (Herbert Creek bordering Shoalwater and Corio Bays Area Ramsar Wetland) exhibited drier conditions than the other sections (Figure 3 and Figure S10 in the Supplementary Materials). A shorter season is apparent with POS occurring in May to July. The NDVI trendline occurred at a lower level than the other sites, but with less temporal variability. Monthly rainfall was 0–390 mm across the time series. Substantial rainfall years in 2010–2011 corresponded with the La Niño pattern of precipitation throughout Queensland but with less intensity (Figure 10). A large decrease in NDVI occurred in 2007 after prolonged low rainfall conditions. The NDVI trendline at this site reflects the El Niño year of 2006/2007 with a low rainfall level and the La Niña year of 2007/2008 with a high rainfall level. POS consistent NDVI years over the 3 sections of our study region were in 2010 with POS in May/June, 2015 with POS in June, and 2019 with POS in June/July.

**Figure 10.** Interpolated, filtered NDVI and rainfall trendlines for the southern section of our study region, Herbert Creek bordering Shoalwater and Corio Bays Area Ramsar Wetland, rainfall data gathered from St. Lawrence Station (station number: 33065).

#### **4. Discussion**

#### *4.1. Change Dynamics in Mangrove Forests*

With one of the richest flora in the world comprising thirty-nine species, old-world tropical mangroves found in the Indo-Pacific including Australia are notable for their attributes of species diversity, abundance, and succession and they are therefore considered to be the most dominant and important mangroves globally [78,79]. We examined decadal changes in mangrove forest cover in Great Barrier Reef catchments using a machine learning approach with the online cloud-computing platform, GEE. With GEE, we compiled a dense time series of NDVI-based Landsat imagery, inputted the time series to the TIMESAT program and investigated the phenodynamics of mangroves at different locations in our study region. This work has relevance to the maintenance of biodiversity, conservation of estuarine wetlands and mangroves, and the role of coastal forests in climate change. To our knowledge, this is the first study to examine and compare the phenology of tropical mangroves with remote sensing time series across different sites adjacent to the Great Barrier Reef, on the east coast of Australia.

We demonstrated that the mangrove extent derived from the RF classification in our study area decreased by 1480 hectares (−2.31%) over the three-year decadal window, 2009–2019. Our findings corroborate a previous study of mangrove deforestation and degradation that described a total of 11,285 hectares lost worldwide from episodic disturbance (climatic events such as dieback, disturbance mortality, climatic extremes, sea level variability) between 2000–2016 with 8004 hectares from Australia [80]. Global annual mangrove loss rates through deforestation, determined from remote sensing datasets, are estimated to be between 0.26% and 0.66% [81]. Decoupling human-induced mangrove mortality such as degradation, from declines due to episodic disturbance events is inherently difficult to achieve [82]. Moreover, spatial inconsistencies can be present as tropical coastal zones are affected by clouds and atmospheric contamination due to their proximity to oceans and the climate regime [83]. Remote sensing-based time series approaches for change detection require a series of images captured continuously throughout a time period. Thus, there is a need for substantially extra and regular image acquisitions over the area of interest [84]. As with previous studies on mangrove change dynamics [26,85,86], we found that the GEE

platform provided efficient data preparation and processing capability for classification and land cover mapping with spectral, auxiliary, and phenology variables.

VIs, particularly NDVI are often used as variables in research related to regional remote sensing forest assessments and shown to be linked not only to canopy structure and leaf area index (the area of single sides leaves per area of soil) but also to canopy photosynthesis [29]. The contrast between strong near-infrared reflectance and low red reflectance of green vegetation is the basis for most VIs. NDVI, and related VIs are functional variants of the simple SR index [87]. We considered four VIs as predictor variables in our classification, SR, NDVI, MNDWI, and GCVI. VIs are commonly used in land cover classification as they increase the classification accuracy, for example: [15] selected NDVI and [88] chose SR to quantify mangrove changes over the Western Arabian Gulf and United Arab Emirates, respectively; [26] generated MNDWI for mapping mangrove extent in China. A major limitation of NDVI however, is that it saturates when mapping high amounts of biomass and cannot distinguish dense vegetation [87]. Another essential point is that NDVI is sensitive to soil, atmospheric effects, and aerosols [29,59]. Furthermore, in open mangrove canopies tidal stage can produce NDVI variations when there is no change in the vegetation amount or condition. Whereas GCVI has been developed as a proxy for chlorophyll content with Landsat surface reflectance data in crops [89], to our knowledge no previous authors have incorporated GCVI for mapping tropical forests or wetlands.

#### *4.2. Mangrove Phenology*

Our multitemporal analysis based on mean-NDVI Landsat time series (from a representative region in the study area, that is, the mid-section for approximately each month in the study period 1999–2019) revealed a strong seasonality in mangrove forests. Specifically, we found a high growth in winter and a trough of growth in late spring to summer. Such cyclical effect in seasonality has been proposed for mangroves in other tropical regions [31] and can be explained by the fact that local climate properties affect the resource availability of water through seasonal influences of tides, rainfall, and river flows and with variation in the strength of ENSO cycles [90]. Considering all available Landsat 8 imagery we demonstrated that the Fourier-style harmonic regression method provides an explicit depiction of temporal shifts in a region with substantial climate variability. Although rarely met in practice, a stringent requirement of investigating continuous intra-annual dynamics of vegetation phenology, is to acquire error-free remotely sensed data at regular time instants [91], and we successfully achieved this goal. The start, length, peak, and end of season are phenological markers that indicate the photosynthetic activity over the growing season. Specifically, we found SOS (start of the season) to be December, EOS (end of the season) to be September/October, and the LOS (length of the season) to be 10 months, regulated by the precipitation pattern in the study region. We suggest that synchronization of new leaf growth with dry season (winter) shedding, shifts canopy composition toward younger, higher radiometric green-efficient leaves, explaining obvious seasonal increases in ecosystem photosynthesis in our study region [92,93]. For instance, studies conducted on *Xylocarpus* spp. phenology-related litter fall in Australia demonstrate that leaf shedding occurred in the mid to late dry season on the tropical east coast [94] and *Avicennia marina* displayed low leafing and shedding rates throughout the year with a seasonal increase in winter [49]. The phenology of *Rhizophora stylosa* in both Northern and Eastern Australia was documented by [95] who detected lower photosynthetic values during the wet season and higher values during the dry season with Landsat imagery, in accordance with our study. Tropical regions experience variable climate conditions with frequent cloud cover and high temperatures during spring and summer months that may limit the capacity of mangroves to photosynthesize. Consequently, greater NDVI values occur in the winter months under cooler temperatures and more radiance. Conversely, a few studies have reported lower spectral index values during the dry season and higher values during the wet season [96,97].

The mismatch that we uncovered between the gradual contraction in the season length, and the increase in seasonal productivity can be explained by mangrove response to the resource availability of water and perhaps nutrients in tropical, coastal soils. Higher rainfall conditions as occurred from the summer of 2009 following six years of drought, produced an associated increase in photosynthetically active growth and biomass production in mangrove forests. High rainfall increases mangrove growth rates by helping to supply sediments and nutrients thereby stimulating productivity [97]. According to [98] water availability influences several ecological processes including vegetative cover, species composition, and community biomass of tidal wetlands. Freshwater is essential for mangroves to maintain the turgidity of cells and tissues, mechanisms of cell augmentation, and stomata regulation [99]. Consequently, in drought years resource-use strategies become geared towards maximum survival as opposed to rapid biomass gain.

#### TIMESAT Analysis

Our TIMESAT-based approach with the adaptive Savitzky-Golay filter generated smoothed data for each time step, and seasonality parameters for each identified growing season over the time series. To provide a clear and interpretable vegetation signal we chose to implement the Savitzky-Golay filter as it iteratively tightens the search window in order to capture very rapid increase or decrease in the data and follows local variations in the seasonal curve more closely [100]. Our findings that a disproportionately bimodal pattern is evident in the time series e g., late summer and winter peaks, confirm studies in phenologybased litterfall conducted on the east coast of Australia with *Avicennia marina* [101], which postulated that bimodality in leaf production is interrelated with position in the inter-tidal zone and frequent tidal inundation. More importantly, abiotic factors such as hydrology and photoperiod and biotic factors such as the timing of the reproductive cycle, appear to be linked to mangrove vegetative phenology [102].

The seasonal dynamics derived using TIMESAT, with mangrove NDVI time series as input, showed differences in the spatio-temporal patterns of growth dynamics. While this is the case, the three different sections of our study region are analogous in that all sections are tropical estuarine ecosystems with the same mangrove species such as the stilted mangrove *Rhizophora stylosa*, the grey mangrove *Avicennia marina*, the yellow mangrove *Ceriops tagal*, and the river mangrove *Aegiceras corniculatum*. Macroclimatic drivers such as rainfall and temperature regimes are recognized for their ability to limit and influence processes within coastal-estuarine wetlands [103]. The premise that rainfall has the greater importance in affecting biodiversity, biomass, and phenology of tropical mangrove vegetation [49] is supported by our study. For example, our observation that peaks and troughs in the NDVI trendline correspond with the trend of monthly rainfall and are magnified by climate events occurred at all three study sections, albeit with differential affects. Accordingly, we found the northern site that was also the wettest had the highest NDVI values and the southern site that was the driest had the lowest NDVI values. The site with the greatest variability in the trendline was the mid-section (Rocky Dam Creek/Cape Palmerston National Park). The mid-section was also the site with moderate monthly rainfall indicating that mangrove growth variability responds to localised patterns.

We also found that the La Niña event of 2010–2011 and its associated rainfall changes impacted vegetation growth, increasing NDVI values at all section sites. Another essential point is that the Southern Oscillation Index spikes in La Niño years, particularly 2010–2011. Indeed, a study conducted at Corio Bay in the Ramsar wetland (a generally euhaline site) in the El Niño year of 2006–2007 reported salinity reached exceedingly high values (up to 40 parts per thousand) towards the mouth. Conversely, at Corio Bay in the La Niña year of 2007–2008 salinities down to inordinately low values (23 parts per thousand) were recorded, which were associated with flooding of the adjacent Fitzoy River and local catchment run off [43]. Notably, our analysis shows that extreme climate-related salinity values such as these, are reflected in the mangrove NDVI trendline where a trough occurs in 2006–2007 followed by rise in 2007–2008 at our southern section (Herbert Creek bordering Shoalwater

and Corio Bays Area Ramsar Wetland). Our results are particularly concerning for tropical mangrove forests as they highlight contemporary research on the intensification of ENSO in the coming decades with climate change [10,82].

Episodic disturbance events such as tropical cyclones and climatic extremes are purported to be responsible for approximately 80% of global mangrove loss from natural causes since 1990 [82]. In addition to coastal damage from strong winds [104], impacts from tropical cyclones include high rainfall that exacerbates the hydrological and hydrodynamic conditions [105]. Moreover, the anomalies in rainfall induce lower salinity, greater turbidity, as well as higher dissolved nutrients and phytoplankton biomass in the coastal waters when compared with typical conditions [106]. The frequency and intensity of cyclones and storms has increased as a consequence of greater sea-surface temperature, with further escalation predicted under most climate change scenarios [107,108]. Greater cyclone frequency has compounding effects on mangrove forests and may extend or halt mangrove recovery [107]. Our findings show decreases in the NDVI values associated with the occurrence of tropical cyclones in each section of our study region. Tropical Cyclone Yasi impacted all three sites in 2011, Tropical Cyclone Jasmine destroyed the northern section site in 2012, and Tropical Cyclone Debbie devastated the mid-section and southern section sites in 2017, which was the highest precipitation event over the time series.

#### *4.3. Limitations of the Study*

Remote sensing data and tools are fundamental methods for measuring land cover, but there are key drawbacks in the change detection of wetlands such as mangrove forests. The need for assessing accuracy of a map generated from any remotely sensed product has become a universal requirement and an integral part of any classification project. The user community needs to know the accuracy of the classified image data being presented [109] otherwise the classification lacks robust validation which could have serious implications for some users and may lessen their confidence in remote sensing as a source of land cover data. Accuracy assessment (1) allows self-evaluation and to learn from mistakes in the classification process; (2) provides quantitative comparison of the classified maps; (3) ensures greater reliability of the resulting maps/spatial information. Therefore, we presented robust validation methods that detail the user's and producer's accuracies of change with Kappa coefficient. Furthermore, we compared the mangrove forest classified maps with Digital Earth Australia high-resolution imagery (at the continental-scale) for visual refinement and to evaluate the performance of our method [56]. The second key drawback is that there can be a similarity of spectral reflectance in the NIR and SWIR regions between mangrove forest and the mixture of other natural vegetation and water at high tide compared to that obtained at low tide [26,86,110]. The third drawback is that quantification is difficult in long-term studies because inconsistencies and gaps in data collection due to observation error are common [111]. The fourth drawback is that temperate climates use simple spring indicators and measures related to 1 January that are not appropriate for tropical phenology because of the circularity of the data. Thus, tropical phenological approaches necessitate taking account of data circularity, be flexible, quantitative, and attribute confidence to conclusions [112]. Finally, a key drawback is that in similar climate conditions, genotypes of mangrove species spread across a study region may vary enough to exemplify different phenophases. In addition, although mangroves are highly specialized and occupy unique and heterogeneous suites of habitats within Australia, they are not a genetic entity but an ecological one [98,112,113]. Consequently, our study has accessed nearby weather stations to our section sites such that localised and subtle climate variations (known to influence phenophases) in location are recorded, and this has reflected in the timeline signals.

#### *4.4. Implications for the Conservation of Mangrove Forests*

Growing evidence of mangrove's ability to sequester and store carbon has brought attention to their conservation value from international climate change policy and decision makers [114]. Although estuarine monitoring is almost completely absent at Shoalwater and Corio Bay Areas Ramsar Wetland, it is critical to establish a monitoring program to provide baseline assessments of vegetation seasonality in relation to weather patterns and detect responses to extreme events. The Ramsar Convention on Wetlands is an intergovernmental treaty with Parties to the Convention from 171 countries. The purpose of the Ramsar Convention is to promote wetland conservation and wise use to ensure that the benefits of wetlands contribute towards meeting the UN Sustainable Development goals, Aichi Biodiversity Targets, Paris Agreement on Climate Change, and other related international commitments [115]. The United Nations' 2030 Agenda for Sustainable Development together with SDGs require the monitoring and mapping of mangroves through the sustainable development indicator 6.6.1: "change in the extent of water-related ecosystems over time" [11]. Monitoring of mangrove forests is critical to preserve the ecosystem services they provide such as pollution control, detoxification, and regulation of natural hazards. Hydrologic conditions, land use, and management are important drivers of coastal vegetation communities with many land use changes developing over decades. So far, studies on tropical forests and mangrove cover with phenodynamics have been limited by the unavailability of suitable remote sensing imagery. Dense time series of satellite images has allowed a comprehensive understanding of mangrove phenology and its relation to climate factors across different sites. Our results are of relevance for conservation strategies of tropical coastal wetlands in Australia that account for foundation species and their seasonal responses to climate drivers.

#### **5. Conclusions**

In this paper, we demonstrated an operational method to (1) map mangrove forests in Central Queensland, Australia using an NDVI time series, which was built in GEE from Landsat imagery, and (2) examine seasonal and inter-annual changes in mangrove phenology extracted with the TIMESAT program, and linked to natural and extreme climate variability. Surprisingly, estuarine monitoring has been absent in the coastal ecosystems of Shoalwater and Corio Bay Ramsar Wetland. Mangrove forests in this southern section of our study region are bounded by an internationally recognized Ramsar Wetland Treaty and therefore changes in their ecology are important for Australia's reporting commitments to the Ramsar Convention. Policy guidelines of the Convention underline that urgent action is needed to raise awareness of the benefits of wetlands, with greater safeguards for their survival. We propose that a regulated monitoring program be instigated so that information can be gathered on the change dynamics and seasonality of estuarine communities and mangroves under threat from episodic disturbances.

Monitoring is essential, ideally before and particularly after an extreme event, to improve our knowledge of the connections among climate, weather, and biodiversity. Our objective to acquire long-term and spatially extensive data was due to the highly variable nature of biological data, and the need to examine periodical events at different locations. The detection of responses to extreme events is important to improve understanding of the connections between climate, extreme weather events, and biodiversity. Such knowledge is essential to inform conservation management attempts to mitigate the impacts of extreme events with ongoing climate change.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10 .3390/rs13153032/s1, Figure S1: Annual sea surface temperature anomaly for the Great Barrier Reef (1999–2020) from Australian Bureau of Meteorology using NOAA ERSST V5 gridded data, Figure S2: Annual rainfall anomaly for Queensland (1900–2020) from Australian Bureau of Meteorology, Figure S3: Linear regression detrended time series of all available Landsat 8 imagery with NDVI for the study region, Bowen to Shoalwater and Corio Bays Area Ramsar Wetland, Central Queensland, Figure S4: Length of season and long integral (total productivity in the annual cycle) in the time series for our study region, Figure S5: The peak NDVI value (maximum recorded NDVI value in the time series) and short integral (the difference between the function describing the season and the base level from season start to end) in our study region, Figure S6: Amplitude (the range between the

maximum NDVI and the base level or the maximum increase in canopy photosynthetic activity above the baseline) and the short integral (the difference between the function describing the season and the base level from season start to end) in our study region, Figure S7: Season year (1999–2020) with seasonal metrics from the Savitzky-Golay filter in TIMESAT, Ampl. (amplitude), L.deriv. (growth rate). R.deriv. (rate of senescence), Figure S8: Southern Oscillation Index and rainfall trendline for the mid-section of our study region, rainfall data gathered from Koumala Station (station number: 33038), Figure S9: Mangrove forest seasonal trend decomposition in TIMESAT with Landsat NDVI data and dense time-series 1999–2020, northern section of our study region, Repulse Bay, black line indicates the seasonal trend, grey dots display the start and end of the growing season, Central Queensland, Figure S10: Mangrove forest seasonal trend decomposition in TIMESAT with Landsat NDVI data and dense time-series 1999–2020, southern section of our study region, Herbert Creek bordering Shoalwater and Corio Bays Area Ramsar Wetland, black line indicates the seasonal trend, grey dots display the start and end of the growing season, Central Queensland.

**Author Contributions:** S.R.P. and D.A.C. jointly conceived and designed this study. H.P.P. provided input into the contents of the paper, suggestions for analysis and editing of the draft and final versions. D.A.C. processed and analysed the data and led the writing of the paper with all authors making contributions. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

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

#### **References**

