An Analysis of the Early Regeneration of Mangrove Forests using Landsat Time Series in the Matang Mangrove Forest Reserve, Peninsular Malaysia

: Time series of satellite sensor data have been used to quantify mangrove cover changes at regional and global levels. Although mangrove forests have been monitored using remote sensing techniques, the use of time series to quantify the regeneration of these forests still remains limited. In this study, we focus on the Matang Mangrove Forest Reserve (MMFR) located in Peninsular Malaysia, which has been under silvicultural management since 1902 and provided the opportunity to investigate the use of Landsat annual time series (1988–2015) for (i) detecting clear-felling events that take place in the reserve as part of the local management, and (ii) tracing back and quantifying the early regeneration of mangrove forest patches after clear-felling. Clear-felling events were detected for each year using the Normalized Difference Moisture Index (NDMI) derived from single date (cloud-free) or multi-date composites of Landsat sensor data. From this series, we found that the average period for the NDMI to recover to values observed prior to the clear-felling event between 1988 and 2015 was 5.9 ± 2.7 years. The maps created in this study can be used to guide the replantation strategies, the clear-felling planning, and the management and monitoring activities of the MMFR.


Introduction
Mangrove forests are woody plants found in the intertidal zones in tropical and subtropical latitudes [1,2]. As such, mangroves are subjected daily to changes in temperature, exposure to water and salt, and different degrees of hypoxia [3]. Mangroves have adapted to these changing conditions and are resilient ecosystems that can additionally cope with environmental changes and disturbances such as storms, hurricanes, mean sea-level fluctuations and climate variability [3]. In areas where conditions have remained relatively stable, mangroves have thrived and natural cycles of succession have prevailed [4,5].
Mangrove restoration and reforestation has garnered interest in recent years due to the carbon sequestration capacity of these forests [6,7]. Silviculture is one of the oldest reasons, among others, to pursue mangrove restoration projects [4,5]. Mangrove reserves such as the Matang Mangrove Forest Reserve (MMFR) in Malaysia, Bintuni Bay Mangroves in West Papua, Indonesia, and the Sundarbans in Bangladesh are examples of mangrove forests that are under silvicultural management [8][9][10]. In these cases of silvicultural management, restoration and reforestation of mangrove forests, it is essential to understand the conditions that favor mangrove propagule dispersal and establishment, mangrove growth and succession, and carbon sequestration.
Remote sensing has been frequently used to monitor changes in mangrove forest cover at local and global levels, e.g., [11][12][13][14][15][16][17][18][19][20][21][22][23]. With an increasing array and availability of satellite sensor data, there is an opportunity to routinely observe mangrove forests and support silvicultural management, restoration and reforestation projects. Using optical data, focus has often been on the use of spectral vegetation indices, such as Normalised Difference Vegetation Index (NDVI), to detect forest disturbances [24,25]. However, indices that exploit the short wave infrared (SWIR) spectra, such as the Normalised Difference Moisture Index (NDMI), have been found to be more effective than the NDVI in detecting forest disturbances [20,[26][27][28]. The use of the SWIR reflectance can provide information about vegetation moisture content and about the upper canopy structure [29]. Moreover, temporal changes in these vegetation indices after a forest disturbance can also be used to detect forest recovery [20,25,26,30].
Studies on time series analysis have the potential to provide information of the past and current changes of mangrove forests [31][32][33][34][35][36][37]. Time series based on vegetation indices have been used to analyse changes in mangrove cover and phenology [32,33,36]. Nevertheless, the use of these types of time series to quantify the mangrove forest recovery after disturbances is still limited [32,33,35]. In this study, we focused on the MMFR located in Peninsular Malaysia, as this reserve is silviculturally managed and provides a unique opportunity to establish whether (i) time series of Landsat sensor data could be used to detect the clear-felling events that regularly take place in the reserve, and (ii) the NDMI annual time series could be used to trace back and quantify the early recovery of mangrove forest patches after a clear-felling event.

Study Area and Silvicultural Management
The MMFR has been under management since 1902 for charcoal and pole production [38]. The reserve is a riverine mangrove forest comprising 27 mangrove species [8]. The reserve has an area of 40288 ha and is divided in four different administrative zones: Protective (17.4% of the reserve's total forest area), productive (74.8%), restrictive productive (6.8%) and unproductive (1%) (Figure 1) [8]. The productive and restrictive productive zones are under silvicultural management and are the only areas where the wood extraction occurs. They are composed of forests dominated primarily by Rhizophora apiculata Blume and R. mucronata Lamk. The silvicultural management consists of a 30-year rotation cycle with two thinnings at 15 and 20 years after clear-felling [8]. The protective zones are composed of diverse mangrove genera such as Avicennia, Sonneratia, Bruguiera and Rhizophora. Patches of dryland forest exist within the protective zones. The reserve provides ecosystem services, including provision of wood for charcoal and pole production, coastal protection, conservation of flora and fauna, and mangrove propagule production [8]. The unproductive zones are represented by lakes and infrastructure areas, including villages, charcoal kilns and administrative buildings [8]. Each zone has a different species composition. The productive and restrictive productive zones are mainly composed of R. apiculata and R. mucronata species. The protective zones are more diverse in terms of tree species, including Avicennia, Sonneratia, Bruguiera and Rhizophora mangrove forests, and dryland forests. For this study, field data were collected in the north of the MMFR (indicated in circles). The grey areas are outside the management of the reserve. Maps adapted from Weidmann et al. [39] and Ariffin and Mustafa [8].
The reserve is located in the west coast of Peninsular Malaysia (Figure 1). The reserve has a tropical climate, with an average air temperature ranging from 22 • C at night to 33 • C during the day [8]. Relative humidity ranges from 99% in the early morning to 60% around noon [40]. Peninsular Malaysia has monsoon season between November and March, and between May and September [41]. The annual rainfall ranges from 2000 mm to 2800 mm [8]. The area is subjected to semidiurnal tides and the tidal amplitude range is 3.3 m [42]. Freshwater flows into the reserve via streams and ground run-off [42].
Every ten years the local management defines the management plan for the reserve [8]. During the period of analysis of this study four different management plans have been developed: 1980 to 1989, 1990 to 1999, 2000 to 2009, and 2010 to 2019 [8,40]. It is noteworthy that the silvicultural practice of a 30-year rotation cycle with two thinnings at 15 and 20 years after clear-felling has remained constant for each plan. Within each management plan, areas to be clear-felled or thinned are pre-assigned and are further divided into smaller areas, called coupes (areas between 2.2 and 6.6 ha), which are then allocated to approved charcoal contractors [8]. Clear-felling activities are performed every year by different contractors, but a coupe will only be clear-felled once every 30 years [8].

Methods
We followed the workflow outlined in Figure 2 to detect the clear-felling events in the MMFR and quantify the early regeneration of mangroves within the logging coupes. Specifically, Landsat data for the period 1988 to 2015 were used to generate a single cloud-free image or a multi-date composite for each year (28 in total). Based on these data, the single Normalised Difference Moisture Index (NDMI) and Normalised Difference Vegetation Index (NDVI) layers were derived. The NDMI time series was found to be optimal to map the extent of clear-felling events and the early regeneration of mangrove forests within coupes. The resulting extent and regeneration maps (i.e., clear-felling and recovery time) were validated against local management plans and a dataset of 135 reference time series. The following sections provide a summary of the processing undertaken.
Remote Sens. 2019, 11, x FOR PEER REVIEW 4 of 18 which are then allocated to approved charcoal contractors [8]. Clear-felling activities are performed 112 every year by different contractors, but a coupe will only be clear-felled once every 30 years [8].

115
We followed the workflow outlined in Figure 2 (Table 1). Surface reflectance Level 2 data products were downloaded from the United States

132
Geological Survey (USGS) Earth Explorer site [43,44]. The surface reflectance data products from 133 TM and ETM+ were generated using the Landsat Ecosystem Disturbance Adaptive Processing

Optical Sensor
Year

Time Series Creation
For each year from 1988 to 2015, an annual surface reflectance image was generated from individual cloud-free scenes or composites of the several scenes within one year. Only one image was used for each year, given that logging coupes are typically cleared over a two-year period and only once every 30 years. To create a multi-date composite image for each year with prevalent cloud cover, we first selected the image with the least amount of clouds and cloud shadow. Where the land surface was obscured, we replaced affected areas with the pixel values from another image of the same year. On average, we only had to replace values in 3 % of data contained in the image with the least amount of cloud cover from 1988 to 2015. The final image was comprised of data from several dates but this was considered acceptable given that clear-felling takes place over an extended period and is not seasonally dependent. The final data stack consisted of 28 images, with each band representing the surface reflectance of a single year. Two additional image stacks representing the NDVI and NDMI (Equations (1) and (2)) for each year were also generated from 1988 to 2015 based on the surface reflectance [29,47]: where NIR is the Near-Infrared band, R is the red band and SWIR is the short-wave infrared band.
In the event of a missing value for a pixel due to cloud cover, the NDVI and NDMI values for the corresponding pixel were interpolated based on the pixel's NDVI and NDMI values from the previous and following year. We only encounter this problem for 2012, where 70% of the image was missing due to a failure in the Scan Line Corrector (SLC) of Landsat 7 [48].

Fieldwork
We acquired Unmanned Aerial Vehicle (UAV) imagery in different forests stands between June and July 2016. We collected data on four different types of forests: A stand that was clear-felled between the end of 2015 and beginning 2016, and stands of 7, 14 and 30 years in age. UAV data were acquired using a DJI Phantom 3 Professional quadcopter UAV and the in-built true color camera, and used to create four orthomosaics of 1ha for each type of forest stand [49].

Time Series Analysis
Using the NDVI annual time series, we created a water mask for the MMFR, which includes all areas of open water in the ocean, rivers, creeks and non-vegetated inundated mudflats. If the NDVI value was zero or negative at any point in the time series, the corresponding pixel was classified as water throughout as this approach accounted for different water levels. Areas classified as water were excluded from further analysis. Additionally, areas of non-vegetated terrestrial land covers (i.e., urban infrastructure and dry soils) were associated with pixels with a negative NDMI value and were also excluded. Finally, we excluded land areas located outside the reserve based on the management zones map ( Figure 1).
Within the remaining area, we detected the clear-felling events using the NDMI time series. We evaluated 210785 time series of pixels from 1988 to 2015. For each series of pixels, we evaluated the following cases:

1.
We detected clear-felling events based on two conditions: a.
We first defined a reference year as the year when the NDMI was below 0.288. This threshold was determined through comparison of the NDMI series of pixels associated with clear-felled and non-clear-felled areas. Areas with dense vegetation (including pre-logged and mature mangroves) exhibited NDMI values around 0.5, which contrasted with areas that were clear-felled which exhibited NDMI values below 0.288. b.
The difference between the NDMI of the reference year and the following year was at least 0.275. We defined this second threshold to guarantee that the drop in the NDMI value was sufficient to correspond to a true clear-felling event. This second threshold was also determined by comparing different series of pixels associated with clear-felled and non-clear-felled areas. The approach followed on from previous studies [30,[50][51][52] that have also used thresholds to analyse time series.

2.
We determined the following values for each clear-felling event: a.
The year of clear-felling b.
The year of recovery. This value was determined as the year when the NDMI value returned to the state prior to clearing [26]. This previous state was defined as the median value of all the points in the series before the clear-felling event minus one standard deviation to account for normal fluctuations in the vegetation index. c.
The recovery time. This time is defined as the number of years that the mangrove forest took to regenerate that is, the difference between the year of recovery and the year of the clear-felling occurrence. If this number was one, we considered this event as noise as mangrove regeneration is not possible in a single year. d.
The drop in the NDMI value, calculated as the difference between the NDMI value before clear-felling and the lowest NDMI value in the time series.
It was not possible to determine the recovery time for clear-felling events that had not restored to pre-clearance levels by 2015. Such areas were typically cleared after 2009 and the forests did not appear to have had enough time to regenerate. Additionally, we did not calculate the recovery time in areas that were already clear-felled in 1988 as there was no knowledge as to whether clearing had occurred in or prior to that year. We

2.
The management zone maps and the logging plans included in the management plans from 2000 to 2009, and 2010 to 2019 [8,40]. First, we compared the existing local management zone map against the results of our clear-felling map, on the assumption that clearing only occurs in the productive and restrictive productive areas (i.e., where wood extraction is officially approved). Second, we compared the clear-felling year calculated in this study against the logging plans outlined in those management plans. These logging plans contain the year when the coupes should be clear-felled.

Time Series Creation
Water areas were more easily discriminated using the NDVI time series as compared to the NDMI time series. Water areas exhibited a negative NDVI value as compared to vegetated areas that had a positive NDVI value (Figure 3b,d). The NDMI provided a measure of the moisture content in both vegetated areas and water (Figure 3a,c) and therefore exhibited similar positive values for both. The NDVI annual time series was therefore used to discriminate water bodies from vegetated areas (see supplementary data Figure S1).   (Figure 3a, 3b), and for forested areas (Figures 3c and 3d). The blue line indicates the 237 mean value for each year and the grey area the standard deviation.

238
The rapid decrease in the NDMI following clear-felling is highlighted in Figure 4a The rapid decrease in the NDMI following clear-felling is highlighted in Figure 4a, with values considerably lower than forests that were not clear-felled between 1988 and 2015 (Figure 4c). The NDMI declines primarily because of the reduction in vegetation amount (and hence structural complexity) and exposure of the underlying surface. The lowest value in the NDMI time series was attained once all vegetation had been completely removed from a clear-felled area, which typically occurred within two years following the clear-felling event [8]. Differences in the NDMI between pre-and post-clear-felling were greater compared to the NDVI (Figure 4b), therefore providing better detection of clear-felled areas. Additionally, we used the NDMI time series to detect urban areas in order to exclude them from our analysis (see supplementary data Figure

Reference to Field and UAV Data
The relationship between the changes in the NDMI time series and the forest structure observed on the field is illustrated in Figure 5. We observed a very closed canopy from seven-year-old to 30-year-old forests ( Figure 5 UAV top view). The contrast with a recently clear-felled area is noticeable, where the remaining trunks and dried branches are still present in the area ( Figure 5, Example 2). Moreover, the first signs of mangrove re-establishment are observed on the ground in the clear-felled areas, with some saplings already present ( Figure 5, ground view Example 2). A view from the ground also confirms the high tree density in young forests of seven and 14-year-old stands ( Figure 5 bottom images, Examples 3 and 4). Although older forests generally have a lower density of trees, the canopy is regarded as closed when viewed from above ( Figure 5, Example 1). It is noteworthy that the forest structure prior to clearing is not the same after the NDMI reaches a value that is similar to that of the pre-cleared forests. Therefore, changes captured by the NDMI time series only reflect the early regeneration of the mangrove forests, as these forests will continue their growth further.

Time series Analysis
In Figure 6, two examples show the year of clear-felling and the year of recovery detected for the NDMI time series of two pixels.
The year of clear-felling for each Landsat pixel contained within the MMFR is given in Figure 7.

273
seven-year-old and (4) 14-year-old stands. Following felling, the soil is exposed and trunks and 274 debris remain. A very dense canopy is then formed after ~7 years of growth despite decreases in 275 density as forests age or are thinned.  (4) 14-year-old stands. Following felling, the soil is exposed and trunks and debris remain. A very dense canopy is then formed after~7 years of growth despite decreases in density as forests age or are thinned.

287
The year of clear-felling for each Landsat pixel contained within the MMFR is given in Figure 7. 287 The year of clear-felling for each Landsat pixel contained within the MMFR is given in Figure 7. The mean recovery time (i.e., the period for NDMI values to return to pre-clearance values), taking into account all the clear-felling events that recovered by 2015, was 5.87 ± 2.67 years (Figure 8). The median recovery time (interquartile range) was five years (4-7 years) and the minimum and maximum was two and 25 years respectively. It is noteworthy that recovery times were uneven across coupes clear-felled in the same year (Figure 8b,c).

297
The mean recovery time (i.e., the period for NDMI values to return to pre-clearance values), 298 taking into account all the clear-felling events that recovered by 2015, was 5.87 ± 2.67 years ( Figure   299 8). The median recovery time ( interquartile range) was five years (4 -7 years) and the minimum and 300 maximum was two and 25 years respectively. It is noteworthy that recovery times were uneven 301 across coupes clear-felled in the same year (Figure 8b, 8c).

Validation Time Series Analysis
Comparison with the 135 reference points (Section 2.2.5) indicated that 100% were correctly identified as clear-felled in the corresponding year and 96% determined the same recovery time as the time defined by visual inspection. Incorrect matching of the remaining 4% (five points) was attributed to noisy data (e.g., a sudden increase or decrease in the signal before the occurrence of the clear-felling, see supplementary data Figure S4) and limited information available on the previous state (clearance history). In these cases, the algorithm overestimated the recovery time by one year in two cases, underestimated by two years in another two cases, and underestimated by five years in one case.
In this last case the signal was noisy and defining the recovery time was not straightforward even through visual interpretation (see supplementary data Figure S4).
Reference to the pre-existing forest management maps (Figure 1) established that clear-felling events had taken place in the productive and restrictive productive zones. We also identified changes in three areas that are protective zones according to the management plan but these changes require further investigation to determine the actual cause of change (Supplementary data Figure S5). Furthermore, we compared the map of the year of clear-felling with the management plan map that prescribes the year when the coupes should be clear-felled. The coupes designated by the local management plan emerge in the map of the year of clear-felling created in this study (see Figure 9). However, the clear-felling events detected in this study occurred on average 2.85 ± 1.15 years after the planned year. It is noteworthy that the management plan is a scheme that is defined before the actual clear-felling. Therefore, the satellite imagery precisely indicates when the actual clear-felling took place (see supplementary data Table S1).
Remote Sens. 2019, 11, x FOR PEER REVIEW 12 of 18 clear-felling, see supplementary data Figure S4) and limited information available on the previous 314 state (clearance history). In these cases, the algorithm overestimated the recovery time by one year in 315 two cases, underestimated by two years in another two cases, and underestimated by five years in 316 one case. In this last case the signal was noisy and defining the recovery time was not 317 straightforward even through visual interpretation (see supplementary data Figure S4).

319
Reference to the pre-existing forest management maps ( Figure 1) established that clear-felling 320 events had taken place in the productive and restrictive productive zones. We also identified 321 changes in three areas that are protective zones according to the management plan but these changes  Table S1).

341
and Sader [28] and Wilson and Sader [27] which also found that the NDMI time series are useful for

Time Series Analysis
The NDMI time series could be used to detect the clear-felling events and quantify the recovery time. These results are in accordance with the studies of DeVries et al. [26], Kennedy et al. [20], Jin and Sader [28] and Wilson and Sader [27] which also found that the NDMI time series are useful for detecting forest changes, disturbances and regeneration. The NDMI is indicative of the moisture content and structural complexity of forests and therefore vegetated areas and water have a positive value. Following clear-felling, the NDMI decreases rapidly as vegetation is removed and the soil is exposed. However, the forests regenerate (with assistance from planting) and the NDMI values return to those typical of the pre-cleared forest and remain similar until the next clearing event. Hence, only the first years of regeneration can be quantified as the NDMI saturates in dense vegetation [54][55][56] despite differences in forest structure. Nevertheless, we were able to identify variations in the recovery times between and within the management-defined coupes.
To the best of our knowledge, this is the first study that quantifies early forest recovery times in the MMFR using annual Landsat time series. Previous studies on mangrove regeneration have focused on quantifying new mangrove cover [32,33]. Although knowledge of new areas colonized by mangrove is certainly important and required to improve the management of mangrove areas, measurements other than area are required to quantify mangrove regeneration in the MMFR. Aziz et al. [57] used nine Landsat satellite images from 1978 to 2014 to analyse mangrove regeneration in the MMFR by classifying the forests as young or mature. Nevertheless, they did not use annual time series nor did they quantify recovery times.
This study generated maps of clear-felling that mirrored the coupes defined by the local management, albeit with some minor differences in the shapes of coupes. This corroborates the findings of this study and opens the possibility for using the methods introduced here to support local management of the MMFR. Clear-felling events were detected on a pixel-by-pixel basis and hence, if one pixel was partially clear-felled and some vegetation still remained within the pixel, it was not possible to detect the clear-felling event as there was an insufficient decrease in the NDMI. This case can occur as the management suggests that about eighteen mature Rhizophora trees should be retained per hectare in the clear-felled areas to serve as a propagule source for recolonization [8]. However, this limitation does not hinder the usefulness of our maps to support the local management, as the area of each coupe is practically completely identified in this study. Overall, we confirm the results of Aziz et al. [57,58] highlighting the usefulness of the Landsat archive to support the local management of the MMFR.
We also identified changes in three protective zones (see supplementary data Figure S5), one of which was also studied by Aziz et al. [57]. We found possible changes in the forest structure captured by a sudden drop in the NDMI time series in areas where Avicennia-Sonneratia, Rhizophora and dryland forests are located. Whilst not in the productive zones, these changes highlight the vulnerability of the managed forest to events or processes that can have implications in the ecosystem services that the protective zones provide for the MMFR [4,57].
Mangrove recovery times were shown to vary from two to 25 years after clearing and also differed between and within areas that were clear-felled in the same year. Mangrove propagule dispersal and establishment are complex processes that are influenced by factors such as currents, water depth, wind, salinity, geomorphology, propagule predation, propagule morphology, light and nutrients availability [59][60][61][62][63][64]. Therefore, further studies are needed to explain the variations in the recovery time that take into account the aforementioned factors and increase our understanding of the drivers of mangrove regeneration.

Implication for the Local Management
The results of this study are new and useful for the local management of the MMFR. The management in the MMFR implements a policy of active replanting. Two years after a clear-felling event, the management assesses the cleared area and if the natural regeneration is less than 90%, Rhizophora seedlings are planted where needed [8]. This policy has an impact on the vegetation cover two years after clear-felling and therefore on the NDMI. The recovery times identified in this study extend beyond the two-year assessment period. Hence, this study can inform the local management on where and when to replant. This study can also identify where past planting might not have been as successful as anticipated and therefore can better inform on future management pathways.
Additionally, the results of this study can also be used to support monitoring activities, redefine current and new administrative forest zones, and guide the planning of clear-felling events.
The identification of changes in the NDMI can also offer a view of changes that occur in both productive zones and protective zones. This information can be used to update current management maps and monitor the forests in the zones where wood extraction does not occur but nevertheless are relevant for maintaining ecosystems services in the reserve, such as propagule source, coastal protection and protection of flora and fauna [8]. Moreover, the availability of free satellite imagery opens the possibility to increase the frequency of the monitoring as mangrove areas are difficult to access [49].
The methodology used in this study can also be used in other mangrove areas that are under silvicultural management such as Bintuni Bay in West Papua, Indonesia [9] and the Sundarbans [10]. Bintuni Bay has a silvicultural management focused on Rhizophora species with a rotation cycle of 25 years and a similar replantation strategy as the one used in MMFR [9]. Therefore, the methodology proposed in this study could also be used by the management of this mangrove area in West Papua. Furthermore, the mangrove areas that are not under silvicultural management could benefit from using the NDMI time series to detect disturbances in forest structure. As recently outlined by Curnick et al. [65] even the smallest patches of mangroves threatened by anthropogenic disturbance serve important ecosystem functions, goods and services and may benefit from early detection when disturbed. As we show in this study, these types of time series have the potential to detect disturbances and regrowth in already existing mangrove areas. Although a modification of the current methodology would be necessary to consider the drivers of other types of disturbances, the potential to track regrowth after a disturbance is useful for mangrove conservation and management worldwide.

Conclusions
This study confirmed that clear-felling events and early regeneration of mangrove forests in the MMFR could be quantified by means of Landsat time series. The NDMI was most useful for detecting clear-felling events that regularly take place in the reserve as part of the local management strategy. Additionally, the NDMI was useful to calculate the early regeneration of mangrove forests, with an average recovery time of 5.87 ± 2.67 years. The patterns of logging were similar to those outlined in the MMFR management plan. Therefore, the mapping provided by this study provides new insights into the recovery of mangrove forests in the MMFR and are a valuable approach to guide future clear-felling management, replanting strategies and monitoring.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/11/7/774/s1; Figure S1: Water mask; Figure S2: NDMI and NDVI time series for urban areas; Figure S3: Time series analysis: NDMI drop calculation; Figure S4: Validation time series analysis by means of comparison with point dataset; Figure S5: Validation time series analysis by means of comparison with the management zones map; Table S1: