TIME SERIES ANALYSIS FOR MONITORING SEAGRASS HABITAT AND ENVIRONMENT IN BUSUANGA, PHILIPPINES USING GOOGLE EARTH ENGINE

Seagrasses are marine flowering plants which are part of a highly productive coastal ecosystem and play key roles in the coastal processes. Unfortunately, they are declining in area coverage globally, and seagrass losses can be attributed to climate change such as sea-level rise, increase in sea surface temperature, and decrease in salinity, as well as human-related activities. The objective of this research is to assess the historical changes in the seagrass habitat and environment of Busuanga, Philippines using time series data available in the Google Earth Engine (GEE) platform. These include satellite data such as MODIS, Landsat 5, 7, and 8, and SeaWIFS. Reanalysis data such as HYCOM was also utilized in this research. Results from HYCOM data show that there has been a 0.0098 °C increase in the sea surface temperature per decade in Busuanga while MODIS data indicates an increase of 0.0045 °C per decade. Moreover, HYCOM data also shows an overall average of 0.76 mm in sea surface elevation anomaly and a decreasing trend in salinity values at 0.0026 psu per decade. Chlorophyll-a concentration has a minimal increase based on results from MODIS and SeaWIFS. Aside from changes in water parameters, changes in the land also affect seagrasses. Forest loss may cause increased siltation in the coastal ecosystem which can lead to seagrass loss. Based on the results of Landsat satellite image processing, there has been forest cover loss in Busuanga with the highest loss occurring in 2013 when super typhoon Yolanda ravaged the island. Lastly, results from the linear spectral unmixing of 778 Landsat images from 1987-2000 show that the average percent cover of seagrasses in Busuanga were declining through the years.


INTRODUCTION
Seagrasses are unique marine flowering plants that can live entirely submerged in water (Short and Coles, 2001). They serve as a nursery for juvenile fishes and as feeding grounds for various marine animals including endangered species such as sea turtles and dugongs (Koedsin et al., 2016;Roelfsema et al., 2009). They help in improving the visibility and water quality in their area and they act as buffers to reduce currents and erosion (Hemminga and Duarte, 2000), as well as increase oxygen levels in the water column. Furthermore, they have various commercial uses and provide high-value ecosystem services (United Nations Environment Programme, 2020;Watson et al., 1993;Waycott et al., 2009). Unfortunately, seagrasses are declining in coverage globally.
Seagrasses can be found in all coastal areas of the world, except along Antarctic shores, and they thrive in coastal zones that are most heavily influenced by humans (Hemminga and Duarte, 2000). However, seagrasses are one of the most rapidly declining ecosystems on earth, with rates of disappearance at 7% per year (Waycott et al., 2009). In the Southeast Asian region, the largest seagrass extent can be found in the Philippines (Fortes et al., 2018). Since 1980, seagrasses have been declining globally at a rate of 110 square kilometers per year and 29% of their known extent has disappeared since 1879 when they were first recorded (Waycott et al., 2009). However, these values may have been underestimated because the study did not include Southeast Asia, which is the global hotspot of seagrasses, because of lack of data (Sudo and Nakaoka, 2020). Recent analysis of seagrass beds in Southeast Asia show seagrass to be declining at an average rate of 10.9% per year since the late 1990s (Sudo et al, under review). Evidence points to human disturbance as the primary cause of seagrass decline globally (Duarte, 2002;Hemminga and Duarte, 2000;Short and Wyllie-Echeverria, 1996), which is exacerbated by climate change (Harley et al., 2006;Short and Neckles, 1999;Waycott et al., 2009). Studies predict more seagrasses will be lost if proper resource management will not be put into place (Borum et al., 2004;Turner and Schwarz, 2006).
Proper management of seagrasses is necessary for their protection and survival. However, this is hindered by lack of information and resources. A better understanding of the seagrass environment is crucial for their conservation and to prevent further loss. Historical data can provide information about the health and status of seagrasses and their environment. Unfortunately, remote areas are less likely to have historical ground monitoring data of the seagrass habitat and environments available. Recent advancements in technology for environmental monitoring, such as remote sensing, may push the proper management of seagrasses in the right direction.
Remote sensing is a tool suitable for large scale mapping and monitoring of the coastal environment. However, the processing of large datasets for environmental monitoring presents some challenges, for instance, availability of expensive hardware required for computationally intensive and time-consuming data processing. Among the recent innovations in remote sensing which can solve such problems is the development of the Google Earth Engine (GEE) platform.
In 2010, GEE, a cloud-based platform for analyzing geospatial data on a planetary-scale (Gorelick et al., 2017), was launched by Google (Wicker, 2010). It creates a stage for the public to utilize Google's vast capabilities in cloud storage and computing. Several geospatial datasets are available in the GEE public data catalog such as Landsat, Moderate Resolution Imaging Spectroradiometer (MODIS), Sea-viewing Wide Field-of-view Sensor (SeaWIFS), HYbrid Coordinate Ocean Model (HYCOM), etc. Lansdat, MODIS and SeaWIFS are products of satellite systems while HYCOM is a data-assimilative hybrid isopycnal-sigma-pressure coordinate ocean model which is a product of collaborating institutions such as the University of Miami, the Naval Research Laboratory (NRL), and the Los Alamos National Laboratory (LANL) (Chassignet et al., 2007). The GEE platform made these products easily accessible in their cloud storage and developed a system for time-consuming and computationally intensive applications like time series analysis of remotely sensed data to be easily accomplished.
The main objective of this research is to monitor the seagrass habitat and environment in Busuanga, Palawan, Philippines using time series analysis of available data in the GEE platform. This includes analysis of sea surface temperature, sea-level rise, salinity, chlorophyll-a concentration, forest loss, and percent cover change of seagrasses.

Study Area
The study area is in Busuanga Island (Figure 1  The coastal ecosystem in Busuanga is relatively unspoiled (Cadigal et al., 2017) in comparison to other areas in the Philippines which are damaged due to human-induced disturbances, bad tourism and aquaculture practices, direct mechanical damages, release of toxic compounds into the coastal waters, etc. Busuanga has thriving seagrass beds, mangroves, macroalgae, and coral reefs (Calamianes IFRM Plan, n.d.).
However, tourism and urban development are starting to increase in this part of the country which is why it is imperative to monitor the coastal resources for protection and conservation (Calamianes IFRM Plan, n.d.). Unfortunately, historical monitoring data is scarce in Busuanga.

Extraction of Time Series Data in GEE
Remote sensing is a valuable tool for environmental monitoring especially in areas where historical data is limited. The GEE platform takes advantage of Google's massive capabilities in data storage and processing. The extraction and processing of data in GEE involve three major stages: selection, filtering, and plotting of time series data for analysis. In the selection stage, a collection of available images is created which is carried out by specifying the start and end date of specific satellite images covering the study area. A shapefile of the study area boundary was uploaded to GEE and was used in the selection stage. The next step was filtering the image collection for specific bands such as sea surface temperature from HYCOM and MODIS, chlorophyll-a concentration from MODIS and SeaWIFS, and sea surface elevation anomaly and salinity from HYCOM. Furthermore, a 2kilometer buffer from the land, as shown in Figure 1, was created, and was used as a boundary for analysis of the time series data. The last step is plotting and downloading the monthly and/or yearly spatial average values of the various parameters within the 2-km boundary. After downloading, time series plots of the data downloaded were created for visualization and further analysis of the data.
A time series is a set of data obtained from observations collected sequentially over time (Meghea et al., 2012). There are two general objectives of using time series analysis: (1) to understand or model a sequence of random variables at fixed sampling intervals and (2) to predict future values based on the history of the series and possibly the driving factors related to the series (Cryer and Chan, 2008). Time series can be separated into trends and seasonal changes that can be modeled deterministically with mathematical functions of time (Nielsen, 2011). Analysing the trend constitutes the use of empirical methods to quantify and explain variations in a system over a period of time (Chandler and Scott, 2011). Time series can be applied in various applications such as physics, economy, engineering, environmental science, etc. In this research, it was applied to monitor and assess the seagrass habitat and environment using remotely sensed data and global models. Examples of the time series data downloaded from GEE are shown in Figure 2 such as sea surface temperature and chlorophyll-A concentration from MODIS (2002-2020) and sea surface elevation anomaly and salinity from HYCOM . The locally estimated scatterplot smoothing (loess), a smooth curve that locally reduces the residual variance (Bărbulescu, 2016), was also included in the time series plots to visualize the general trend of the data. Furthermore, the time series were decomposed using the additive model to determine the trend and seasonality of the parameters.
Unfortunately, observation data of water quality parameters in Busuanga are not available. To validate the trend of the extracted values, the MODIS sea surface temperature was compared to values from HYCOM while the chlorophyll-a concentration from MODIS was compared to SeaWIFS extracted chlorophyll-a concertation values within overlapping time frames of the available data.

Forest/Land Cover Change Monitoring
GEE made available in their platform the global forest loss data created by Hansen et al. in 2013. In their study, they defined forest cover as 25% or greater canopy closure at the Landsat pixel of 30 meters by 30 meters spatial resolution for trees greater the five meters in height. They focused on the conversion of forest cover to non-forest cover wherein 2000 is the reference year for pairing Landsat images (Hansen et al., 2013). The forest loss data in Busuanga from 2001-2018 was extracted in GEE by getting the spatial average within the boundaries of the area. This information was plotted as time series and downloaded for further analysis.
Furthermore, to confirm the forest loss in the study area, Normalized Difference Vegetation Index (NDVI), one of the most widely used vegetation indices in remote sensing (Campbell and Wynne, 2011), were also extracted from available MODIS data from 2001-2018. Monthly spatial NDVI averages were calculated and then the trendline was determined.

Monitoring of Seagrass
Remote sensing has been extensively utilized through the years for seagrass monitoring using various sensors and platforms. In a study by Traganos et. al. in 2018, they presented the capabilities of the GEE platform for large scale seagrass mapping which can eventually lead to a global-scale monitoring. In their research, they used support vector machine to classify Sentinel-2 satellite images and achieved an overall accuracy of 72% (Traganos et al., 2018). In this study, to monitor the status and changes in the seagrasses of Busuanga, a linear spectral unmixing method was performed using Landsat 5, 7, and 8 images. The Landsat satellite imagery was chosen for this research because of its extensive collection since 1980s which is crucial in determining the possible seagrass loss. Spectral unmixing was the method used to classify seagrasses because there are no available historical seagrass training data for Busuanga which can be used for supervised classification methods. Only field data points from recent field surveys were available and these were utilized to determine the spectral signatures of the sample points. Then, the extracted spectral signatures were used to identify similar spectral signatures from archived images using the spectral umixing method. Spectral unmixing is a sub-pixel classification wherein the proportions of different classes contributing to the spectrum measured at a particular pixel are estimated (Rees, 2013). Five classes were used in the spectral unmixing: seagrass, corals, macroalgae, sand, and deep water. A total of 778 Landsat satellite images of the study area with less than 20% cloud cover from December 1987 to April 2020 were available in GEE and utilized in this research. Among them, 296 were from Landsat 5, 278 from Landsat 7, and 204 from Landsat 8. Shown in Figure 3 is the total number of available images per year. Pure spectra or endmembers were gathered from images within August to December 2019 which corresponds to the field survey dates conducted in Busuanga. Seagrass percent cover data from the field survey were utilized for validation of the seagrass spectral unmixing results. After obtaining satisfactory results from the validation, the endmembers were used for the linear spectral unmixing of all images from 1987-2020. The seagrass field data survey was carried out using the SeagrassWatch (Seagrass-Watch, 2010) manual wherein three 50m transects were laid out parallel to each other and perpendicular to shore, and 25m apart, on the each site, and seagrass percent cover were estimated at 5 meter intervals within 0.25 m 2 quadrats along the transects. The average values were computed for each site to represent the seagrass percent cover for the area. Due to the closeness of some sites, values were averaged from them as well and seven sites (shown in Figure 13) were used to validate the seagrass percent cover results from spectral umixing.

Increase in Sea Surface Temperature
Global mean sea temperature continues to rise year after year, and extreme temperatures will have negative effects on seagrasses (Short and Neckles, 1999). Depending on the species, the effects of increased temperature will vary. Minimal increase in temperature may have positive effects, such as increased productivity, on some species while extreme change may result in variations in seasonal and geographic patterns of species and abundance (Lee et al., 2007).

Sea Level Rise
Based on a study by Saunders et al in 2013, there may be a 17% loss of seagrasses globally if the sea level rises by 1.1m in 2100 (Saunders et al., 2013). An increase in water depth results in attenuation of light intensity penetrating the water which causes lower seagrass productivity (Short and Neckles, 1999). However, this may be compensated by a landward migration of seagrasses (Brodie and De Ramon N'yeurt, 2018;Short et al., 2016). The problem arises when there are infrastructures in the coastal zone blocking the migration of seagrasses which will lead to their mortality. Reanalysis data from HYCOM ( Figure 6)

Salinity Change
Climate change causes more intense storms and increased rainfall which leads to increased freshwater run-off from the land to the coastal ecosystem (Short et al., 2016) causing a decrease in salinity. On the other hand, sea-level rise causes greater saltwater intrusion in the coastal and estuarine locations which may result in an increase of salinity values. Salinity change, whether an increase or decrease, will impact seagrasses (Short and Neckles, 1999). Different species have varying responses to salinity changes. Effects include a decrease in productivity and biomass, limitation in reproduction and distribution, changes in community structure, or even mortality (Duarte, 2002;Short and Neckles, 1999). Based on HYCOM data from 1993 to 2020, a decreasing trend of 0.0026 psu per decade in salinity values occurs in Busuanga, as shown in Figure 7, which also indicates seasonality with the highest salinity happening in May while the lowest salinity occurs during August to October. This coincides with the dry and wet season, correspondingly, in Busuanga. During the wet season, higher precipitation occurs which causes a decrease in sea water salinity.

Chlorophyll-a Concentration
Chlorophyll-a concentration is useful in estimating productivity and it is also an indicator of light stress in seagrasses (Gallegos et al., 2004). An increase in chlorophyll-a concentration may indicate the presence of an increased phytoplankton population in the water which may result in less light penetrating the water and reaching the seagrass beds. Fortunately, results from MODIS (2002( -2020( ), and SeaWIFS (1997, as shown in Figure 8, indicate a minimal increase in chlorophyll-a concentration in Busuanga with only 0.00013 mg/m 3 and 0.000037 mg/m 3 per decade, respectively.

Forest Loss
Forest loss will result in increased siltation in coastal ecosystems (Duarte, 2002;Milliman and Meade, 1983). Increased siltation will result in increased light attenuation and possible burial of seagrasses which may lead to a decline in seagrass diversity, biomass, and production (Duarte, 2002;Terrados et al., 1999).  As indicated in Figure 10, the highest loss occurred in 2013 when super typhoon Yolanda ravaged the island. It is also significant to note that in 2006, an airport was established in Busuanga for cargo planes and small aircrafts. It was renovated in 2008 to be able to accommodate larger passenger planes. Based on the graph, a large increase in forest loss happened after the establishment of the airport because of growth in development and tourism on the island. Figure 11. Map of forest loss in Busuanga (2000Busuanga ( -2018 Looking at the map of forest loss (Figure 11), there was no concentrated area of loss. Patches of forest loss can be seen dispersed within the study area which indicates no concentrated massive development. Aside from using forest loss data, vegetation indices such as Normalized Difference Vegetation Index is a good gauge for land cover change. Based on the MODIS NDVI time series in Busuanga from 2001 to 2019, as shown in Figure 12, there is a decreasing general trend in NDVI values. This trend confirms the loss of forest cover in the study

Seagrass Percent Cover Change
Seagrass area is in decline all over the world because of natural and anthropogenic factors. It may be difficult to determine the exact cause of the decline because it is usually a combination of various factors. In this study, different factors which may lead to seagrass loss were presented. To determine if there is a decline in seagrass cover in the study area, the linear spectral unmixing method was performed on archived Landsat images in GEE. A map of the result of the linear spectral unmixing of seagrasses using Landsat 8 images from August to December 2019 is shown in Figure 13. Overlaid on this map are the location and percent cover of seagrasses gathered from the field around Busuanga.

Figure 13. Map of seagrass percent cover in Busuanga
To validate the results of the linear spectral unmixing of seagrass percent cover, a correlation analysis was performed using the field survey data. Figure 14 shows the results of the validation indicating an R square of 0.7741. Figure 14. Validation of seagrass spectral unmixing using seagrass percent cover gathered from field surveys After validation, the same method was applied to a total of 778 Landsat 5, 7 and 8 images from 1987 to 2020 to determine if there is seagrass loss in the study area. All of the images used were calibrated top-of-atmosphere (TOA) reflectance available in GEE. To summarize the spectral unmixing results, a 500-meter buffer was created from the land and was used as a boundary of pixels and all values within this boundary were averaged.
The plots of the yearly average of the linear spectral unmixing results of seagrasses is shown in Figure 15 and it indicates a decreasing trend in seagrass percent cover in Busuanga. The high seagrass percent cover in 1992 may be due to a significantly low number of available images for this year which may have resulted to amplification of errors. Water column correction was not implemented in the data processing because the parameters needed to accomplish this vary from image to image such as depth and attenuation coefficient. However, because the values are averaged annually, this lessens the errors in the time series like seasonal, tide and water quality variability when compared to processing single images per year. Figure 15. Yearly average seagrass percent cover from spectral unmixing

Correlation Analysis
The annual average of seagrass percent cover, sea surface temperature, seas surface elevation anomaly, salinity, chlorophyll-a concentration, and forest loss were calculated, and a correlation analysis was performed between these parameters. The correlation results are shown in Figure 16. Seagrass percent cover was positively correlated with sea surface temperature, salinity, chlorophyll-a concentration from MODIS and forest loss while it was negatively correlated with sea surface elevation anomaly and chlorophyll-a concentration from SeaWIFS. The opposing correlation between seagrass percent cover and chlorophyll-a concentration from MODIS and SeaWIFS is due to the difference in the range of available data between them. SeaWIFS data was available from 1997-2010 while MODIS data was available from 2002-2020. It is also important to note that there was a minimal increase of chlorophyll-a concentration in Busuanga through the years. The correlation between seagrass percent cover and the parameters were relatively low because seagrass loss can be attributed to a combination of various phenomena. The highest correlation with seagrass percent cover is sea surface elevation anomaly which means that as the sea surface elevation anomaly increases, the seagrass percent cover decreases based on the data presented. Furthermore, relatively high positive correlation was observed between sea surface temperature and forest loss, sea surface temperature (HYCOM) and chlorophyll-a (MODIS), and salinity and forest loss. The high correlation does not prove causation, but it shows that there is a relationship between the variables. Positive correlation indicates that as one variable increases, the other also increases based on the data shown in this study. Figure 16. Correlation matrix between the annual average of seagrass percent cover, sea surface temperature (SST_HYCOM and SST_MODIS), sea surface elevation anomaly (SSEA), salinity, chlorophyll-a concentration (Chla_SEAWIFS and Chla_MODIS) and forest loss

CONCLUSIONS
The seagrass habitat is part of the coastal environment which is among the most productive ecosystems on earth (Duarte et al., 2018). It serves as shelter to a variety of marine species and provides valuable goods and services to humans (Koedsin et al., 2016;Roelfsema et al., 2009;United Nations Environment Programme, 2020). The health of this habitat is highly influenced by coastal water quality. Degraded water conditions leads to loss of seagrass habitat (Waycott et al., 2009). For that reason, it is important to monitor various coastal water quality parameters, as well as land cover changes, for the benefit of the marine organisms and people living in coastal communities that rely on seagrasses.
GEE provides a platform that makes available several geospatial datasets suitable for large scale mapping and monitoring of the coastal environment. It also eliminates the need for expensive hardware and large storage for processing such datasets. Furthermore, GEE has made it efficient to process timeconsuming and computationally intensive applications, such as time series analysis of remotely sensed data.
In this study, using time series analysis of datasets in Busuanga using GEE, we found out that sea surface temperatures were increasing at 0.098°C per decade based on HYCOM data, while MODIS data showed an increase of 0.0045°C per decade. Sea surface elevation anomaly also increased at a rate of 0.0027 mm per decade, while salinity was decreasing at 0.0026 psu per decade, based on HYCOM data. On the other hand, chlorophylla concentration increased minimally through time, at a rate of 0.00013 mg/m 3 and 0.000037 mg/m 3 based on MODIS and SeaWIFS data, respectively. Using the Hansen et. al. global forest loss data available in GEE, we observed forest loss in Busuanga. Lastly, based on spectral unmixing results of 778 Landsat images from 1987-2000, we observed seagrass percent cover decline in Busuanga. The correlation analysis showed that seagrass percent cover was positively correlated with sea surface temperature, salinity, chlorophyll-a concentration from MODIS and forest loss, while it was negatively correlated with sea surface elevation anomaly and chlorophyll-a concentration from SeaWIFS.
For future research, we suggest employing advanced statistical methods to further analyze the relationship between water quality parameters and seagrass loss in Busuanga. Also, to improve the accuracy of seagrass monitoring using satellite images, water column corrections may be applied. However, this may be a challenge because of the extensive collection of images needed to be processed and the parameters needed for water column correction will differ from one image to another.