*3.1. Data*

#### 3.1.1. Reference Data

Very high-resolution reference data were used for the recent epoch. This dataset is available as a MAXAR Vivid basemap within the ArcGIS software [67,68]. These cover the study area with data from November 2009 to January 2017. About 90% of the study area is covered with 46-cm-pixel data from GeoEye-1 (10 December 2010, 16 December 2011, 3 January 2013, 17 December 2013, 10 April 2014, 8 January 2015), 60-cm-pixel data from QuickBird-2 (11 February 2010, 3 October 2010, 12 June 2013), and 50-cm-pixel data from WorldView-2 (1 December 2011, 16 February 2013, 13 January 2014, 12 March 2015, 17 December 2015). Thanks to the familiarity with the study area, the broad land cover classes that were targeted in this paper were relatively easily identifiable on the very high-resolution imagery. This was also the case for the degraded mangroves, which presented the additional advantage of being spatially confined within the coastal zone, in general, and the mangrove system, in particular.

#### 3.1.2. Landsat Data

The choice of Landsat data was driven by the need to coincide with as many other NDR studies as possible, so that comparisons could be drawn between them. Two such studies were identified: the one by Ayanlade and Drake [23] and the study by Kuenzer et al. [27]. The latter was particularly targeted, as it is the only one that has attempted to map the "degraded mangrove" class. The choice of the three epochs was also driven by the availability of the reference data and the SAR imagery.

We used all the dry season (December to February) Level 1 surface reflectance Landsat 4, 5, 7, and 8 images centred around 1988 (±2 years), 2000 (±2 years), and 2013 (±2 years) with less than 80% cloud cover from the USGS EROS Data Center for the eight WRS-2 tiles covering the study area (path 187, row57; p188, r55; p188, r56; p188, r57; p189, r55; p189, r56; p189, r57; p190, r56). Only the non-thermal bands were used, and clouds and cloud shadows were removed using F-mask [69,70]. Finally, the Normalised Di fference Vegetation Index (NDVI) [71] was calculated. From the resulting 7-band image stacks (i.e., six non-thermal bands, plus the NDVI), spectral-temporal variability metrics were calculated [33,34,72,73]. For the recent epoch, five statistics for each of the seven bands were calculated: the standard deviation, the mean, and 3 percentiles (25th, 50th, and 75th). This brought the total layers for this epoch to 35. However, as data availability for the first two epochs was problematic (Figure 3), we limited the number of statistics per band to 2 (mean and st. dev.) and the total number of layers to 14.

**Figure 3.** Number of available observations from the Landsat USGS Level 1 archive for (**a**) the first epoch; (**b**) the middle epoch, and (**c**) the more recent epoch.

#### 3.1.3. Radar Data

Radar data were chosen for testing whether their addition to the optical metrics could improve the land cover classification. For the recent epoch, we employed the 2015 global 25 m resolution L-band Synthetic Aperture Radar data from the Advanced Land Observing Satellite-2 (ALOS-2) PALSAR-2 sensor via Google Earth Engine's API. The data are free and open access with two polarisations (HH and HV) and are currently available for 2015 to 2018. To increase the utility of the SAR data, we used Google Earth Engine to calculate a series of Gray-Level Co-Occurrence Matrix (GLCM) texture variables [72]. GLCMs are a series of localised texture metrics that quantify the statistical properties of a layer over a moving window [74]. We calculated seven GLCM layers (mean, variance, homogeneity, contrast, dissimilarity, entropy, and second moment) [75]. These statistics were calculated over both 3 × 3 and 9 × 9 windows, resulting in 15 layers per SAR backscatter (one backscatter + seven 3 × 3 GLCM layers + seven 9 × 9 GLCM layers), totalling 30 layers for the year 2015.

For the middle epoch, we acquired JAXA's 25 m resolution JERS-1 tropical region mosaics for the year 1996, the only year that such data are available over the Niger Delta Region. One polarisation is available (HH), from which we calculated 15 GLCM layers to use in the classification.

#### *3.2. Land Cover Mapping*

#### 3.2.1. Sampling and Validation

In total, 185,504 samples were taken for the epoch centred around 2013. For the first and second epochs (i.e., 1988 and 2000), TimeSync-Plus v4.6 was used [76] to check for unchanged pixels at the 2013 sample locations. This resulted in 142,045 and 99,220 samples, respectively, for which we could confidently say that no change in the Landsat time series occurred. During classification, half of these samples were used for training and half for validation.

#### 3.2.2. Image Classification & Post-Classification Processing

We developed the land cover classification using Random Forests classification models. Random Forests have been used successfully to classify Landsat imagery, thanks to their effective handling of correlated predictors and reduced tendency toward overfitting [77]. We used the 'RStoolbox' and 'randomForest' packages within the R statistical environment [78]. One optical only model was tested for the first epoch, while for the middle and most recent ones, we tested the performance of optical only and optical + SAR metrics (Figure 2). Based on the accuracies achieved, the outputs from the best performing models were chosen for the middle and more recent epochs. A 3 × 3 majority filter was applied to the outputs from all 3 epochs to ge<sup>t</sup> rid of the 'salt and pepper' effect of the classification. Finally, based on our knowledge of the study area, expert rules were applied to correct for some classification errors [72].

## *3.3. Intensity Analysis*

Aldwaik and Pontius [52] devised a methodology that characterises patterns of land change quantitatively. It provides a mathematical framework that compares a uniform intensity to observed intensities of temporal changes among land cover classes (or 'categories') [79]. There are three levels of analysis, with each level exposing different types of information given the previous level of analysis. The first level, i.e., the interval level, examines how the size and speed of change vary across time intervals. The intensity of the rate of annual change is estimated using the following equations [52] (for notation, see Table S1 in the Supplementary Material):

$$S\_t = \frac{area\,\,of\,\,change\,\,during\,\,interval\,\,[\,\,Yt\,\,,\,\,Yt+1\,]\,\,area\,\,of\,\,study\,\,region}{\,\,duration\,\,of\,\,interval\,\,[\,\,\,Yt+1\,]} \times 100\%\,\tag{1}$$

$$
\Delta U = \frac{\text{area of change during all intervals/area of study region}}{\text{duration of all intervals}} \times 100\%. \tag{2}
$$

The second level is called "category level" and it examines how the size and intensity of gross losses and gross gains in each land cover class vary across classes for each time interval. This level identifies which land cover classes are relatively dormant or active in each time interval [52]. Equations (3) and (4) provide the intensity of a class' annual gain and loss, respectively:

$$\mathbf{G}\_{lj} = \frac{\text{area of gross gain of class j during } [\text{Yt}, \ \text{Yt} + 1] / \text{duration of } [\text{Yt}, \ \text{Yt} + 1]}{\text{area of class j at time } \text{Yt} + 1} \times 100\%,\tag{3}$$

$$L\_{\rm ti} = \frac{\text{area of gross loss of class I during [Yt. } \text{Yt+1]} / \text{duration of [Yt. } \text{Yt+1]}}{\text{area of class i at time } \text{Yt}} \times 100\%. \tag{4}$$

The third level, the "transition level", examines how the size and intensity of land cover class' transitions vary across the other classes that are available for that transition [52]. At each level, the method tests for stationarity of patterns across time intervals and identifies which land cover transitions are particularly intensive in a given period. Aldwaik and Pontius [52] provide a detailed explanation of the limitations concerning where the transition from a particular land cover class *m* to a class *n* can occur. For example, if a given land cover class *n* exists at a particular location at the initial time, then class *n* cannot gain at that place. If class *n* gains, then it must gain from locations that, initially, are not class *n*. If class *n* gains uniformly across the study area, then this class will gain from other classes, in proportion to the initial sizes of these land cover classes. Alternatively, class *n* might intensively avoid gaining from some particular class(es) and might intensively target gaining from some other class(es). Given the observed gross gain of class *n*, Equations (5) and (6) identify which other classes are intensively avoided versus targeted for gaining by class *n* in a given time interval:

$$R\_{\rm int} = \frac{\text{area of transition from i to n during } [\text{Yt. } \text{Yt} + 1] / \text{duration of } [\text{Yt. } \text{Yt} + 1]}{\text{area of class i at time } \text{Yt}} \times 100\% \tag{5}$$

$$\mathcal{N}\_{\text{tn}} = \frac{\text{area of gross gain of class } n \text{ during } [\text{Yt }, \text{ Yt} + 1] / \text{duration of } [\text{Yt }, \text{ Yt} + 1]}{\text{area that is not class n at time } \text{Yt}} \times 100\%. \tag{6}$$

We used the *intensity.analysis* package in R to carry out the processing (https://cran.r-project.org/ web/packages/intensity.analysis/vignettes/README.html).

#### *3.4. Landscape Pattern Analysis*

Post-classification comparison is most informative about changes in the composition of a landscape but gives us little—only visual—information about the spatial characteristics of these changes and the distribution of landscape elements. Landscape pattern analysis using landscape metrics provide us with additional information about the structure of changes, such as landscape fragmentation and patch aggregation or dispersion, as well as their changes in time. With the latter, we can observe changes in landscape spatial configuration through time.

We followed the approach used by Gounaridis et al. [53] and selected a number of class-level metrics [80] in order to study the changes in the spatial configuration and patterns of the 'mangrove' and 'degraded mangrove' land cover classes. We used 'Percentage of Landscape' (PLAND) as a measure of class abundance, and the 'number of patches' (NP), 'landscape patch index' (LPI), and 'patch area median' (AREA\_MD) to study fragmentation of the classes of interest. With regard to patch shape analysis, we used the 'area weighted mean patch shape index' (SHAPE\_AM), and for the aggregation of these classes, we used the 'area weighted mean Euclidean nearest neighbour distance index' (ENN\_AM) along with its standard deviation (ENN\_SD). Finally, we also used the aggregation index of 'percentage of like adjacencies' (PLADJ). Table 1 provides a listing of the selection of landscape metrics used in this study, together with a short description of their correlation with mangrove forest fragmentation. For more information, refer to McGarigal and Marks [80] who provide a full description of the metrics, including their mathematical formulas.


**Table 1.** Selection of landscape metrics used in this study with a short description of their relationship with mangrove forest fragmentation.
