**4. Discussion**

Extracting phenological metrics of mangrove forests from satellite images is an ongoing field of research. We contribute to this field by (1) presenting a novel, data-driven method to extracting phenology from satellite imagery; and (2) applying the method in evergreen forests across Australia. We also expand this field of research by presenting the dual phenology of mangrove forests, as described in the literature. By demonstrating that there is more than one period of leaf growth in these ecosystems, we also highlight the need for deeper investigation into the detection and causes of dual phenology in evergreen forests and the need for long term field studies that validate satellite observations. Importantly, with the use of satellite images, we have demonstrated that plots with similar species can have different phenologies. Phenology is, in turn, site-dependent and hence should not be described using a single logistic or sinusoidal curve.

Because phenology is site-dependent, methods like the one presented here can be used with long term imagery archives to determine the baseline of mangrove phenology in stable ecosystems, and compare it to stressed ecosystems. Early detection of ecosystem stress, indicated by changes in phenology, could lead us to detect, predict, and hopefully prevent, mangrove dieback events.

#### *4.1. The Phenology of Rhizophora Stylosa*

Authors have described, in situ, the phenological traits of *R. stylosa* around the world, however, few have attempted to compare mangrove phenology across regions. [16], for example, used satellite

imagery and a sinusoidal model to describe seasonal variations of mangrove forests in the Yucatan peninsula in Mexico. They found that the spectral index values are lower during the dry season and higher during the wet season, an outcome that di ffers from our findings. Across all our study sites, we found lower EVI values during the wet season and higher EVI values during the dry season. These di fferences may be due to the geographical location and the species composition of both studies. *R. mangle*, *Laguncularia racemosa*, *Avicennia germinans* and *Conocarpus erectus* dominated their study site, while *R. stylosa, A. marina and C. tagal* dominated ours. Despite *R. stylosa* being the dominant species in our site, several species of mangroves may contribute to the apparent phenology of each pixel given the resolution of the Landsat images (i.e., 30 m). Determining the contribution of each species to the apparent phenology is an ongoing field of research.

Another big di fference between our study and that of [16] is that the species that dominate mangrove forests in their study sites show a unimodal phenology response across the Yucatan peninsula. In contrast, we found bimodal phenology signals in our field study site, as well as in the sites described in the peer-reviewed literature (Figure 4). This bimodal phenology is not new [29], and later [6] described *R. stylosa* and other mangrove species as having two distinct periods of leaf growth each year. The novelty of our analysis has allowed us to demonstrate that these two periods of leaf growth can be detected from satellite imagery using semi-parametric methods. While not all mangrove ecosystems display dual phenology, GAMs create the capability to detect it, giving researchers a tool for inductive reasoning and opening the doors for further mangrove research to validate satellite observations.

Moreover, the method presented here accounts for di fferences in leaf growth intensity between the sites (e.g., Figure 6C,E). Resolving these two periods of leaf growth is important to improve current models of carbon and water fluxes and their relationship to climate forcing. Descriptions at the plot level certainly help us understand the bimodal response of EVI over time, however, the number of leaves that a tree produces might not be the only explanation for a bimodal response in EVI.

Recently, [46] suggested that the seasonal response of EVI in an evergreen forest was bimodal due to layers of the canopy responding in opposing seasonal patterns. When there is a decrease in leaf area index of the upper canopy, the lower canopy takes advantage of the extra light to increase its leaf area. Following the findings by [46,47], propose that the seasonal response of EVI comes from variations in leaf area index and photosynthetic capacity (i.e., younger, more e fficient leaves), rather than from climate and weather patterns alone. Indeed, leaf ontogeny, demography and longevity could influence satellite observations, with leaves of di fferent ages having varying amounts of chlorophyll, carotenoids, and water, resulting in slightly di fferent spectral signatures. These assertions need to be tested in mangrove ecosystems, to determine if the bimodal response of EVI is related to the canopy structure or net leaf production.

In this study, we have demonstrated that seasonal and inter-annual changes in leaf gain and net leaf production are related to changes in EVI. [32] suggested that leaf production of *R. stylosa* in Northern Australia is most evident during the wet season (December through May), however, this species produces new leaves throughout the year. Similarly, [33] found that *R. stylosa* has the highest values of leaf production and leaf fall between December and April. The authors also indicated that this species has a net leaf gain between January and August, and a net leaf loss between September and December, which coincides with upward and downward trends in the apparent phenology for that site (Figure 6E). Despite the strength of the apparent phenology-net leaf production relationship, when using satellite imagery, other factors a ffect the phenology response of mangrove forests.

Environmental and biological factors such as cyclones, rainfall and tree age are known to alter the phenology and spectral response of mangrove ecosystems [16,48,49]. With the help of satellite imagery and GAMs, we can now look at inter-annual changes in mangrove phenology (e.g., Figure 4) and relate them to environmental or biological factors. Because these factors may change one by one, or several at a time, inter-annual predictions of mangrove phenology are di fficult. A way to improve the prediction of seasonal and inter-annual phenology is to include past and present observations of these factors during the creation of the GAMs. GAMs create a numerical relationship between each factor and the phenology model to determine how influential is the former over the latter. This is an ongoing avenue of research.

Similarly, biological variables such as species composition, growth rate and forest maturity may also affect the spectral response of mangroves, and thus, any model derived from satellite sensors. A mature, dense forest will have a different response to a newly-planted mangrove patch or one that is recovering from a natural or manmade disaster [50]. Because relating 30 years of satellite observations to biological processes requires that in situ data is collected frequently and over long periods, the need for long-term (five or more years) monitoring sites in mangrove forests is evident. Long-term monitoring of mangrove ecosystems is important, especially when relating variations in spectral indices, to changes in tree growth and temperature such as those presented by [51,52].

#### *4.2. GAMs vs. Parametric Methods*

Our ability to detect and forecast mangrove phenology improves our understanding of the ecosystem [1], and our approach greatly differs from others, more commonly used [53]. For example, [54] derived a mathematical function that resembles the phenological phases of forests in the northeastern United States. They aimed to use satellite images to monitor vegetation dynamics at the landscape level, and one of their biggest achievements was that the method did not require any fixed constants or thresholds to be applicable. This method has been used at the local [55] and global scales [23], however the premise that a parametrized mathematical curve fits every plot of land has remained unchanged. The method derived by [54] assumes that the selected model of phenology is: (1) correct, (2) already known and (3) is invariant through space and time [39,56]. Even recent studies [17] insist on these assumptions when applying smoothers and filter to the data prior to detecting the phenology without considering that the data they discard may provide insights into the phenomenon they are trying to model (i.e., phenology). This in itself is not a limitation of the parametric models, but of the analysis workflow selected by the authors. We, on the other hand, used semiparametric GAMs as an estimation method to describe the changing relationship between mangrove phenology and EVI. GAMs let the data define the shape of the phenology curve, allowing for bimodality and skewness to be detected and modelled [56].

As shown above, the relationship between EVI and mangrove phenology changes with space and time. The main limitation of parametric approaches is that those methods are constrained to the particular models evaluated (e.g., [54,57,58]). That is to say, parametric methods assume that the shape of the phenological curve remains invariant and only variations in the frequency and amplitude of the signals are allowed. As demonstrated by [46,47] leaf demography and ontogeny can change the way we detect phenology from remotely sensed data, and as such, we need models that can tell us more than just seasonal variations in greenness. GAMs can fill this gap.

In contrast to parametric methods, GAMs use the data itself (in this case EVI values from satellite imagery) to determine the shape of the relationship. Because the relationship between the predictor and response variable is unknown beforehand, GAMs apply a series of smooth functions and use the data to determine which function is the best fit for a given dataset. This data-driven approach enabled us to demonstrate three key things: (1) the phenology response of mangroves forests dominated by *R. stylosa* is not unimodal, but often bimodal across our study sites; (2) the second leaf growth phase varies in intensity depending on the site. For example, the second leaf flush is much lower in New South Wales (Figure 6E) than in Queensland and the Northern Territory (Figure 6B,D). The reasons behind these differences are not fully understood but may be due, in part, to differences in air temperature, rainfall and water temperature; (3) the phenology response of mangrove ecosystems is site-dependent and GAMs allow us to see small seasonal and inter-annual variations that are otherwise impossible to detect using fully parametric methods. However, to provide a better, more accurate, description of mangrove phenology we need ecophysiological descriptions that span more than 18–24 months to compare to satellite-derived phenology.

#### *4.3. Validation of the GAMs*

Peer-reviewed literature of phenological models often lacks a description of the methods used to validate the models, making it di fficult to compare the performance of the GAMs versus other approaches. Despite this limitation, we validated our apparent phenology using three di fferent, independent sources of information and found that GAMs are good tools to extract the phenology response of mangrove forests from satellite images.

Linear regressions between Observed and apparent phenology values showed a good model fit, with R<sup>2</sup> > 0.40 in half of our study sites despite variations in the raw EVI values (Figure S1, see supplementary materials). Correlation values between the GAM and in situ (Table 3, and Figure S1) data were low (as expected), but there are valid reasons for this. The study by [26] focused on three sites in the Gladstone region in Queensland and aimed at examining potential bioremediation strategies in case an oil spill hit the Queensland coast. The authors never intended to use the field data to validate satellite imagery, hence the di fficulties in correlating one with the other. Furthermore, the dataset consisted of one value per date per site i.e., only three data points per date, and some dates had no values (e.g., April 1997), which reduced the potential for correlation even more. Lastly, the study only gathered data over an 18-month period, limiting our information to one and a half growing seasons. Having only one full season of information limits the number of links we could create between the apparent phenology and the field data. We need longer field studies to understand fully the phenology curves extracted from satellite imagery and to examine the response of mangrove forests to changes in weather and climate patterns.

The validation of the apparent phenology with peer-reviewed literature provided us with two important pieces of information: (1) our models correlate well with the leaf gain intensity and net leaf production reported in the literature, regardless of the year in which those data were acquired; (2) apparent phenology has a two- to three-month lag with leaf gain intensity in most of our study sites. The former is important because it highlights the usefulness of GAMs. The latter tells us that, from a biophysical perspective, EVI responds to the canopy elements that absorb red light for photosynthesis and scatter near infra-red light, while field phenology traces leaf formation and drop. The delayed response of EVI is expected for two reasons: (1) the time it takes newly formed leaves to reach their maximum size, and (2) the net leaf production varies throughout the year. In Figure 8 we show an example of this. At "t1", leaves are scarce and bud breaking, net leaf production and chlorophyll content are low but positive and the satellite captures mainly the background of the mangrove tree (i.e., exposed soil and understory water). At "t2", leaves are growing and new leaves are bud breaking. At "t3", net leaf production peaks, meaning that there are many leaves growing and bud breaking, the satellite captures mainly green material, including chlorophyll, thus EVI is high as well. When "t4" arrives EVI is still high, more leaves are dropping than bud breaking, but leaves from "t3" are still growing and reaching maturity (i.e., high chlorophyll content) hence the lag between peak EVI and peak leaf gain/production. Finally, the tree loses more leaves and the background in the satellite images starts to show again ("t5") and the cycle starts again. The implications of this delayed response of EVI need to be explored because the growth/recovery rate after a natural disaster may not be as evident from satellite imagery as previously thought [50].

As demonstrated here, GAMs can be used in di fferent locations, and with di fferent species. Single species mangrove ecosystems are rare, and it is common to have a variety of species, especially within a 30 m pixel. Each study site here (Figure 1) was chosen due to the dominance of *R. stylosa*, however, the associated species varied from study to study. From the Northern Territory to the New South Wales coast, our study covers di fferent regions of the Australian coastline and uses data from di fferent points in time, thereby we have demonstrated that GAMs are easily transferable through location and time. As a result, GAMs can be used in local, continental, and global scale models of phenology spanning years or decades. We contribute to the remote sensing community by demonstrating that GAMs can be used in combination with remotely sensed data and presenting an alternative to fully parametric models.

**Figure 8.** Time difference between peak leaf production and peak EVI during a given year for a simulated mangrove tree.
