**1. Introduction**

Around the world, the e ffects of changing plant phenology are evident in many ways: from earlier and longer growing seasons to altering the relationships between plants and their natural pollinators [1–3]. Remote sensing techniques allow us to detect subtle changes in plant phenology, and here we present a novel approach to describe phenological cycles of mangrove ecosystems. We contribute to a better understanding of mangrove phenology by investigating physical changes of mangrove ecosystems and how the evidence of change is captured by satellite images. Accurately modelling and predicting mangrove phenology will help us understand not only the seasonal variations but also the long-term trends in the natural cycles of these forests. New models, such as the one presented here, will advance our understanding of how drought, heatwaves and other extreme weather events a ffect mangrove health and growth. Similar to using sea temperature to predict coral bleaching events, we could use phenology to predict mangrove dieback events akin to those of 2015 and 2016 in the Gulf of Carpentaria in northern Australia.

Phenology is related to the life cycle events of plants and animals and their relationship to climatic and other abiotic factors [3,4]. Plant phenology also plays an important role in the carbon cycle in the form of sequestration and storage. Phenological cycles of plants ensure that leafing, flowering and fruiting events occur during the most appropriate season to achieve maximum growth and reproductive success. Mangrove phenology is often described at the species level by relating the time of year when trees flower, fruit or defoliate with suspected drivers like temperature and rainfall [5,6]. For example, [7] described the phenology and distribution of *Avicennia marina* mangroves along the Australian coastline, and [8] described the flowering and leafing phenologies of mangroves in the Darwin region. While these descriptions provide a very valuable baseline for comparison, they often lack the spatial extent and frequency needed for phenological studies [9].

We can monitor mangrove phenology using remote sensing or we can collect in situ data. The main advantage of in situ monitoring is that it provides information at the tree and species level, where observations can be very detailed over a wide range of variables. However, in situ monitoring is challenging, time-consuming and variation in methods and survey e ffort can make it di fficult to compare results [10]. In contrast, the remote sensing approach provides information at the landscape and continental scales and is consistently acquired over time and space [9]. While many studies used space-borne sensors to map mangroves at the global [11,12], continental [13,14], and local scales [15], few have used these sensors to monitor mangrove phenology. [16] used MODIS (Moderate Resolution Imaging Spectroradiometer) data between 2000–2014 to detect mangrove phenology using four di fferent spectral indices in the Yucatan peninsula in Mexico. Similarly, [17] compared the phenology of mangroves to that of the surrounding forests using MODIS imagery. While the temporal resolution of the MODIS sensors is very high (1–2 days), the spatial resolution is coarse (250–500 m). Landsat satellites o ffer a better spatial resolution at the cost of a lower temporal resolution. Despite this tradeo ff, the Landsat archive is key to using remote sensing to monitor mangrove phenology as it provides more than 30 years of imagery at a spatial resolution of 30 m x 30 m and a temporal resolution of 8–16 days [18].

To date, most studies on plant phenology have used fully parametric models, mainly in the form of double logistic or sinusoidal functions [19–21]. These functions may perform well in deciduous or temperate forests, where there is a single, well defined period of leaf production, and a single, well defined period of leaf senescence [22,23], but these methods may not be well suited for mangroves and other evergreen forests. When detecting phenology, one of the main limitations of parametric models (e.g., logistic functions) is that they fail to detect asymmetric trends in leaf growth or senescence [22]. Considering that the growing season of some evergreen forests consists of two periods of leaf growth and death, fully parametric (or model-driven) models have the potential to oversimplify the phenology of these ecosystems. Other, more complex models have also been used to examine plant phenology, mainly in the form of artificial neural networks, however, these methods are known to have mostly been used in croplands [24]. Semi-parametric (or data-driven) models, on the other hand, may be better suited for this task as they do not necessarily assume that there will be a single peak or trough in leaf growth or death. Rather, semi-parametric models use the data to determine the shape of the phenology.

Studies have documented the dual phenology of mangroves [6] and other evergreen forests [25] in the field, but this dual phenology has not been recorded using satellite imagery. Dual phenology refers to two periods of leaf growth; unlike deciduous forests which grow their leaves during the spring, some evergreen forests have two periods of leaf growth every year. In mangrove forests these events have never been documented using satellite imagery, probably due to the use of fully parametric models to detect phenology, or because in situ data collection focuses mainly on litterfall rather than leaf production. The novelty of this study is that we use a semi-parametric method to model mangrove phenology and in doing that we present, for the first time, these two distinct periods of leaf growth described in the literature.

Generalized Additive Models (GAMs) are commonly used in ecology and climate sciences, to examine non-linear relationships between response and independent variables. Here we present a novel, data-driven method to extracting mangrove phenology from a series of Landsat images. We use discrete observations of mangrove forests (i.e., satellite images) and Generalized Additive Models (GAMs) to create a continuous curve of phenology over time, without assuming a certain shape, amplitude, or frequency. Our aims are to (1) use a semi-parametric approach (GAMs) to examine if seasonal changes in biophysical variables are related to seasonal changes in the spectral reflectance of mangrove forests; (2) compare the satellite-derived phenology with a set of field observations and measurements; (3) compare the satellite-derived phenology to peer-reviewed literature describing the phenology of mangrove forests, and (4) determine how the Enhanced Vegetation Index (EVI) responds to leaf gain, leaf fall or net leaf production in mangrove ecosystems across Australia.

This manuscript is organized in the following way: firstly, we describe the site and methods used to collect the field data. Then we describe how we use the literature to create a proxy for mangrove phenology. Afterwards, we describe the use of GAMs and satellite images to detect the apparent phenology of mangroves across northern Australia. Having done this, we present the models of apparent phenology and compare them with the data collected in the field, and the proxies from the literature. Finally, we discuss the results, limitations, and future work.

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

We selected six study sites across Australia to evaluate mangrove phenology from satellite imagery using GAMs. One site corresponds to field observations collected in the Gladstone region (Queensland) in the late 1990s, and the remaining sites (*n* = 5) correspond to qualitative data extracted from peer-reviewed publications (Figure 1). In this section, we first describe the field site followed by the peer-reviewed studies, the image acquisition process, and the time series analysis using GAMs. Finally, we describe the phenology model validation.

#### *2.1. Field Site Description*

The Gladstone region in central Queensland is home to over 100,000 ha of intertidal wetlands, out of which 30% are mangrove forests [26]. The annual mean temperature ranges between 18.4 ◦C and 27.5 ◦C, the mean annual rainfall is 874 mm, and the semi-diurnal tides often range from 1.5–3.5 m but reach up to 6 m. In this region, mangroves of the *Rhizophora, Avicenna*, and *Ceriops* genera are among the most common [5,27]. The field data were collected between July 1996 and August 1998 in two plots located in Fisherman's Landing and one plot on Curtis Island (Figure 2). The sites were dominated by *Rhizophora stylosa* trees and located within the Landsat World Reference System 2 (WRS-2) path 91 and rows 76 and 77.

**Figure 1.** Workflow and location of the study sites used to validate the phenology model. \* shows the location of the field site [26] and two other published studies used.

#### *2.2. Field Observations and Measurements*

The data were collected by [26] in the following way: in each plot, the authors selected mature *R. stylosa* trees between 4–9 m tall and tagged 21 leafy shoots in the upper two meters of the canopy totaling around 378 shoots. The authors conducted monthly inspections and recorded the shoot growth, number of leaves, reproductive parts and number of branch shoots. To measure the amount of litterfall, the authors suspended litter traps (1 m<sup>2</sup> in area) under the selected trees. The traps were suspended above the high tide mark and litter was collected, sorted and weighed on a monthly basis. During the data collection period, the average tide height was 2.49 m according to the historical record for the site [28]. Importantly, [26] never intended to validate satellite imagery with their data, therefore (1) no spectral data were collected, making these data completely independent from the satellite-derived phenology, and (2) we do not anticipate a high correlation between satellite-derived phenology and the field data. The data consist of the mean monthly values of six phenological variables, however, we selected the three that were more relevant for our study (see ([26], pp 105–108) for details). The selected biophysical variables are detailed below:

**Leaves lost [leaves** × **m**<sup>−</sup>**<sup>2</sup>** × **day**−**1]:** The number of leaves that fell into the litter traps.

**Leaves gained [stipules** × **m**<sup>−</sup>**<sup>2</sup>** × **day**−**1]:** The number of stipules that fell into the litter traps. This variable serves as a proxy for the number of leaves produced in a tree.

**Net leaf production [leaves gained-leaves lost]:** The difference between leaves gained and fallen leaves. This measure is an indication the net balance of leaves in the canopy with more or less leaves as leaves appear or fall, leaving the canopy in either debit (=stressed), credit (=growth) or neutral condition.

**Figure 2.** Location of the field sites and mangrove patches in the Gladstone region, Queensland. Aerial images of the study site for 1996, provided by the State of Queensland (QAP5402131/47).

#### *2.3. Published Literature on the Phenology of R. stylosa*

To compare the apparent phenology (i.e., from the GAMs) to other sources of information, we gathered a set of peer-reviewed papers that included *R. stylosa* as target species. The reasons for selecting this species were twofold: (1) it is common throughout northern Australia; (2) a number of studies have described its phenology over a wide geographic area across the Indo West Pacific region. We looked for papers that had a graphical interpretation of leaf fall and/or leaf gain over time and we found six examples (Table 1). We used the published graphical interpretations of leaf fall and leaf gain to determine, in a qualitative way, the times of the year where most leaves grew or fell. All published graphs show "Time" on the horizontal axis and a measure of leaf fall or gain on the vertical axis. In each study, we divided the vertical axis into five equidistant categories (i.e., very low, low, medium, high, very high) and recoded the category for each month (not shown). Finally, we calculated the net leaf production from each study by subtracting the leaf fall from leaf gain values and then compared the three variables with the apparent phenology of each site (Table 1).

From Table 1 one can see that the studies by [29,30] predate the time where Landsat imagery was collected. Between 1974 and 1989, cyclones Dawn (March 1976), Keith (January 1977), Gordon (January 1979), and Kerry (March 1979) affected Hinchinbrook Island, Gladstone, or Proserpine in Queensland. Neither [29] nor [30] mention the effects of cyclones, drought on their respective study sites, either because the field campaigns happened before the cyclones, or there was no significant damage to the trees. While there is no certainty about the effects of extreme weather events for those studies, here we assume that the phenology observed by the authors did not change until satellite imagery was acquired.


*Remote Sens.* **2020**, *12*,

 4008

#### *2.4. Landsat Image Acquisition and Processing*

Digital Earth Australia holds a copy of the Landsat archive (1987–present) for the whole of Australia [18,34]. Digital Earth Australia provided all L1T images (*n* = 668) of our sites between 1987 and 2006 for the Landsat 5 (TM) and Landsat 7 (ETM+) sensors, and performed (1) geometric, (2) atmospheric and (3) Nadir-adjusted Bidirectional reflectance distribution function Reflectance (NBAR) corrections following [35]. All corrections were performed within Digital Earth Australia, using the python programming language. Digital Earth Australia uses the 'Pixel Quality Assessment' algorithm [18] to remove pixels with clouds and cloud shadows, as well as all missing pixels resulting from the Landsat 7 Scan Line Corrector failure; these pixels were removed from the datasets and not used (i.e., there was no gap-filling for ETM+ data). The main reason for not filling the missing pixels is that we wanted to use the GAM with only true values, rather than including additional uncertainty to the GAMs. Having done this, we calculated the Enhanced Vegetation Index (EVI) [36] for each pixel in each image. We chose EVI because studies show that this spectral index does not saturate with high vegetation densities [36], it is better suited than other indices for discriminating vegetation fraction in mangrove ecosystems [37], and it is commonly used for phenology investigations [16,17].

Our study leverages the high temporal density of the Landsat archive, and the overlapping footprints of two or more Landsat scenes, hence increasing the number of usable pixels in a given area. On average, the di fference between the dates when field data was collected and the closest satellite image acquired is 5 days, as shown in Table 2. To compare the GAMs with peer-reviewed literature, we selected a period equal to the time of the data collection plus and minus one year, thereby ensuring that the models had enough input data. In the cases where studies were dated before 1987, we used the first three years of available imagery of the area to create the GAM (Table 1). We estimated the location of the studies from the site descriptions in each publication and created a region of interest of approximately 17 ha of mangrove forests surrounding the study area. Afterwards, we selected only the pixels that corresponded to mangrove forests using the "Mangrove Canopy Cover" product developed by Geoscience Australia ([13]) and applied the GAM to every pixel within our region of interest. This approach ensured that we captured the phenology of the mangrove community instead of a small plot.


**Table 2.** Comparison between field data collection dates and satellite image dates.

#### *2.5. Time Series Analysis Using Generalized Additive Models*

With all images pre-processed, georeferenced and sorted by time of acquisition, we proceeded to create a model of phenology for every available pixel in our study sites using GAMs. Contrary to linear additive models, GAMs are statistical models in which the relationship of predictor and response variables is captured by smooth functions instead of coe fficients [38]. Equations (1) and (2) show the respective linear and generalized additive relationships between one response variable (*Y*) and two predictor variables (*Xi*) for *i* observations [39]:

$$Y = \beta \alpha + \beta\_1 X\_{i1} + \beta\_2 X\_{i2} + \varepsilon \tag{1}$$

$$Y = \beta\_0 + f\_1(X\_{i1}) + f\_2(X\_{i2}) + \varepsilon \tag{2}$$

Noticeably, there is no change in the form of the model, however, there is no assumption that the relationship between predictor and response variables is linear. In Equation (1) an additive linear relationship between *Y* and *Xi* is captured by the slope terms β1 and β2, while in Equation (2) the additive relationship is captured by the "smooth" functions *f*1(·) and *f*2(·). The shape of the "smooth" functions ( *fn*(·), also known as "splines" or "smooth splines", is determined during the computation in an iterative way and can take many forms (e.g., polynomial, linear, quadratic) [39,40]. Each 'smooth' function ( *fn*(·)), or spline, is comprised of several basis functions (*bn*), their coe fficients (β*n*), and where *K* determines the maximum complexity of each smooth function:

$$f(\mathbf{x}) = \sum\_{k=1}^{K} \beta\_n b\_n(\mathbf{x}) \tag{3}$$

As explained by [41], GAMs use several types of splines to determine the relationship between each predictor variable ( *Xi*) and the response ( *Y*) is evaluated at every data point in an iterative way, for every predictor variable (see [42] a detailed description). In this case, and following [43], the base function is a Fourier basis with 2N parameters δ = [*<sup>a</sup>*1, *b*1, ... , *aN*, *bN*] *T*, that allow the construction of a matrix of seasonality vectors for each past and future value of *t*. Importantly, we are neither fitting a Fourier series to the data, nor assuming that the relationship of EVI and time is symmetric. Here, the apparent phenology of each pixel results from the Fourier basis expansion evaluated at each data point, and adding the weighted basis functions. In summary, the detecting apparent phenology is a curve-fitting exercise rather than a time series decomposition one.

Another characteristic of GAMs is that measurements do not need to be evenly spaced in time [43]. This works well in our case for two reasons: (1) pixels with clouds, shadows and other errors are flagged as invalid observations leading to time series with random gaps in both length and timing; and (2) areas where the footprints of two or more scenes overlap will have more observations than areas with no overlap.

To detect mangrove phenology from a satellite-derived data series, we used the Python programming language and the "Prophet" package (version 0.3.post2) developed by Facebook [43]. Facebook designed this package to analyze user engagemen<sup>t</sup> with the social network at di fferent time scales, and to investigate how periodic events such as holidays a ffect that engagement. Similar to mangrove phenology, user engagemen<sup>t</sup> on the social media platform is a ffected by regular and irregular events such as weekends (i.e., regular events) and public holidays, which change every year [43]. In a similar fashion, mangroves are a ffected by regular changes in temperature and rainfall (i.e., seasons) and irregular events such as cyclones or drought. While the time scales may di ffer, the concept of tracing an event (e.g., phenology or user engagement) over time remains the same. We selected this package due to its ease of use and re-purposed it to extract seasonal variations in greenness and phenological metrics from satellite images of mangrove forests.
