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

: 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 out-weighing those of dryland ecosystems per hectare. However, land use change and water regulation infrastructure have reduced connectivity in many river systems and with ﬂoodplain and estuarine wetlands. Mangrove forests are critical communities for carbon uptake and storage, pollution con-trol and detoxiﬁcation, and regulation of natural hazards. Although the clearing of mangroves in Australia is strictly regulated, Great Barrier Reef catchments have suffered landscape modiﬁcations 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 reﬂectance products, accessed in Google Earth Engine. Two land cover classes were generated (mangrove versus non-mangrove) in a Random Forest classiﬁcation. Mangroves decreased by 1480 hectares ( − 2.31%) from 2009 to 2019. The overall classiﬁcation accuracies and Kappa coefﬁcient 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 ﬁndings 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.


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 CO 2 eq by the end of the century [4]. For reference, emissions from the global transportation sector are projected to rise to 11,900 Tg CO 2 eq by 2100 [5]. Although the rate of mangrove loss has decreased substantially since the 1990s, from~2% to <0.4% per 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][9][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][21][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][24][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][27][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.

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 •

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 com-pared 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.

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 (ϕ): where p t is the value of the harmonised NDVI within a pixel obtained from the satellite data at time t, β 0 is a constant component of the regression that designates the start of greenness at each pixel, β 1 is a linear coefficient (slope/linear trend), β 2 is the cosine term coefficient, β 3 is the sine term coefficient, e t 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][71][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 R 2 ). 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/).

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.

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].

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).

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).

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

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.  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 2 = 0.8182), and season (year) and season length (Pearson correlation coefficient r value = −0.91, Adjusted R 2 = 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). Table 3. Savitsky-Golay filter with phenometrics extracted from TIMESAT for the Landsat dense time series of mangrove forest in the mid-section of our study region, Startt.: start of season is the beginning of measurable photosynthesis in the mangrove canopy (month of the year); Endt.: end of season is the end of measurable photosynthesis in the mangrove canopy (month of the year); Length: duration of photosynthetic activity (days); Baseval.: the average of the minimum values left and right of the season (NDVI); Peakt.: time of maximum photosynthesis in the canopy (days); Peakval.: the maximum recorded NDVI value (NDVI); Ampl.: the maximum increase in canopy photosynthetic activity above the baseline (NDVI); L.deriv.: growth rate (day/NDVI); R.deriv.: rate of senescence (day/NDVI); L.integral: total productivity (day/NDVI); S.integral: seasonal productivity (day/NDVI); Startval.: value at start of season (NDVI); Endval.: value at end of season (NDVI), Bowen to Shoalwater and Corio Bays Area Ramsar Wetland, Central Queensland.

Seas.
Startt   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 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. 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

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 TIME-SAT 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

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.

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.

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.

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.

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.