Coral reef applications of Landsat-8: geomorphic zonation and benthic habitat mapping of Xisha Islands, China

ABSTRACT

Being one of the most significant and valuable coral reef systems in the South China Sea, the Xisha Islands has undergone rapid transformation due to increasing stressors from human impacts and climate change in recent years. However, as indispensable information for coral reef monitoring and management, the detailed reef extent, geomorphic zonation, or benthic composition of the Xisha Islands is not well documented. Considering limited access to the Xisha Islands, the rapid development of optical remote sensing technology provides us with a feasible mean for coral reef observation. This study adopted a water depth substitution index – probabilistic inundation (PI) – combined with depth-invariant index (DII) to achieve reef extent exploration, geomorphologic and benthic habitat types classification with unsupervised classification algorithms based on Landsat-8 time-series satellite data. Compared with two open-access datasets, the extent of each independent reef extracted from PI exhibited higher similarity with the actual boundary conditions displayed in RGB (Red-Green-Blue) composite images from Landsat-8. Based on PI and derived slope, we obtained geomorphic zonation classification results, and similarly benthic compositions were retrieved based on PI, DII, and reflectance. The overall accuracy of geomorphic zonation and benthic habitat classification results were 72% and 86%, respectively. We also interestingly discovered that corals of the Xisha Islands may be capable of an ability to resist chronic heat stress as a growth trend of reef area after two successive stress events in 2014–2015 were observed at most reefs. The proposed mapping framework of this study provides a repeatable and flexible scheme in depicting the comprehensive situation of coral reefs at Xisha Islands based only on publicly available remote sensing data without complicated pre-set parameters, which could be easily extended to coral reef research around the world. Simultaneously, the findings also provide requisite information supporting the sustainable management and conservation of coral reef ecosystems in the Xisha Islands.

KEYWORDS:

• Coral reef
• Landsat
• time series
• classification
• Xisha Islands

1 Introduction

In recent years, enhanced optical remote sensing approaches and analytical capabilities have substantially overturned the way we perceive our planet, enabling human management of planetary resources from a large temporal and spatial perspective (Runting et al. 2020). Yet for coral reefs, these improving potentialities to collect and analyze data are typically not synchronously transferred into effective information for environmental preservation for several reasons (Hedley et al. 2018; Hughes et al. 2018; Madin, Darling, and Hardt 2019). Concretely speaking, limited access to field data collection and varied water covering situation apparently restrict mapping accuracy. On the other hand, transformation of remotely acquired observations of coral reefs into user-friendly spatial information for further monitoring, planning, and management demands a consistent and repeatable processing scheme to discretize continuous data regarding natural phenomena into spatial units of manageable information (Kennedy et al. 2021).

High spatial resolution remote sensing imagery (<10 m), coincident field data, and object- or pixel-based mapping methods have been successfully applied to obtain coral geomorphic zonationand benthic composition maps on some typical coral reefs (Andréfouët et al. 2006; Hedley et al. 2018; Pang et al. 2021; Phinn, Roelfsema, and Mumby 2012; Roelfsema et al. 2013; Rowlands et al. 2012; Tian, Zhu, and Han 2020; Xu et al. 2016; Yang et al. 2016). Nevertheless, the application of high-resolution polar-orbiting spaceborne images is often constrained by its exorbitant expense, low temporal resolution, and limited observation duration, which is not conducive for long-term, large-scale and systematic repeated coral monitoring. In contrast and as a compromise, moderate resolution imagery (10–30 m), such as Landsat series with regular repeat sampling, is free for public, providing extensive cost-effective worldwide coverage, which can be traced back to the 1980s (Hedley et al. 2016). The Millennium Coral Reef Mapping Project (MCRMP) successfully depicted ~6000 km² of reefs across Caribbean-Atlantic, Pacific, Indo-Pacific, and Red Sea, proving the potentiality of Landsat series in coral reef detecting and habitat mapping (Andréfouët et al. 2006). Compiled from multiple data sources, including the MCRMP, United Nations Environment Programme World Conservation Monitoring Centre (UNEP-WCMC) published the global dataset of warm-water coral reefs distribution, and it mainly relies on Landsat-7 ETM+ (UNEP-WCMC, WorldFish Centre, WRI, and TNC 2021). In addition, a number of independent studies achieved coral reef investigation with Landsat series imagery from regional to local scale. Roelfsema et al. (2002) have used Landsat TM imagery to investigate the spatial distribution of microalgae on reefs in 1999. Vanderstraete et al. (2006) examined changes in coral reef compositions (coral, sand, seagrass, and macro-algae) from 1987 to 2000 based on a Landsat-5 TM image and a Landsat-7 ETM+ image in the coastal zone of the Red Sea. At the same location, El-Askary et al. (2014) continued this work with a Landsat-8 image added, detecting obvious macro-algae increasing and coral decreasing. Likewise, Palandro et al. (2003, 2008) demonstrated shallow-water coral habitat alteration in the Florida Keys National Marine Sanctuary using Landsat-5/TM and 7/ETM+ images. Hedley et al. (2018) also stated that Landsat-8-based retrieval of bathymetry (R² = 0.8) and benthic distribution (overall accuracy of 49% using six categories) showed acceptable accuracy, though more detailed features could be captured with higher spatial resolution sensors.

Despite the existence of numerous studies on mapping coral geomorphic zonation or benthic compositions, these proposed processing schemes were closely related to specific coral reefs and surrounding water environment. Accuracy and credibility of their results, to some extent, depend on remote sensing data and in-situ measured data acquired, which in turn confines the migration application for long-term coral reef monitoring demand in designated areas. In these circumstances, a prompt systematic and universal strategy to conduct long-term and large-scaled surveys for diverse coral reefs is still lacking. This is an especially urgent need for monitoring coral reefs in a timely manner around those remote rarely inhabited islands, such as the Xisha Islands. Yet, regarding rare accessibility to the remote islands and few public in-situ data, a deficiency or even absence of field data in the spatio-temporal domain may lead to the dilemma of investigation scheme validation. Based on such consideration, aiming to provide an instrument to reef restoration, protected area management, and disaster risk reduction, the Allen Coral Atlas has been built by a dedicated team of scientists, technologists, and conservationists (Atlas 2022). It offers a series of comprehensive and unprecedented maps (e.g. benthic habitat maps, geomorphic zone maps) of global coral reefs derived from high-resolution satellite imagery. Benefitted from Allen Coral Atlas, the gap between earth observation data and geo-ecological knowledge of reefs is bridged, as well as the trade-off between applicability at global scales and relevance and accuracy at local scales (Kennedy et al. 2021). Furthermore, it facilitates the following bathymetry/geomorphic/benthic mapping exploration in results’ verification and evaluation.

The Xisha Islands own the largest land area among the four offshore archipelagoes in the South China Sea (SCS), and are home to abundance coral reefs, providing fisheries, coastal protection, and important ecological values (Zuo et al. 2017). Therefore, the health level of the Xisha Islands has become a major national and international concern. Due to a combination of coral disease, crown-of-thorns starfish outbreaks, changes in water temperature and overfishing, not only the coral cover showed a sharp decrease by 21%, but also the coral species composition and community forms of coral species in the Xisha Islands changed obviously over the past 20 years (Li et al. 2018; Pang et al. 2021). In spite of increasing threats and pressures, the Xisha Islands, to some extent, remains one of the least comprehensively mapped coral reef systems in the world. Such deficiency of spatio-temporally explicit coral geomorphic zonation or benthic composition maps further hinders effectively coral monitoring and well-timed coral bleaching warning of the Xisha Islands. Thus, taking a systematical investigation to detecting geomorphic zonation and benthic habitat situation of the Xisha Islands is vital and imperative for disaster risk assessment and marine environmental protection.

Though the Allen Coral Atlas provided relevant coral information of the Xisha Islands, the high-cost of high-resolution satellite imagery, unavailable in-situ validation data, irregular update frequency, and unstable regional accuracy limited its application in dynamic monitoring of regional coral reef conditions. Presently, it should be considered as a reference atlas, rather than precise navigation for long-term reef monitoring. Besides, most researches about coral reefs of the Xisha Islands were focused on several isolated reefs, lacking overall significance (Pang et al. 2021; Wang et al. 2019 ; Wu, Yang, and Ying 2022; Xu et al. 2016; Yang et al. 2016). Therefore, mapping of overall reef extent, geomorphic zonation, and benthic composition of the Xisha Islands based on freely available remote sensing data with highly operable method facilitates for continuous updating remains a task to be completed. The goal of this study is to build a set of complete and simple data processing methods to achieve reef extent, geomorphic zonation, and benthic distribution mapping around Xisha Islands by integrating long-term Landsat-8 images. Furthermore, this study also preliminarily explores the potential applicability of such systematic approach in detecting the multi-year and long-term alteration trend of coral reef area. The results are expected to serve as a practicable and reasonable scientific coral reef mapping framework and technical tool to update the distribution and spatio-temporal variation information of the reefs around the Xisha Islands effectively and efficiently. Eventually, the outcome will realize the accurate translation of outcome obtained based on remote sensing data from producers to end-users (reef managers, policy makers, conservation practitioners), contributing to a future management, development, and protection of the reefs.

2 Materials and methods

2.1 Case study locations

Coral reefs within the SCS are principally distributed in Dongsha Islands (Pratas Islands), Xisha Islands (Paracel Islands), Zhongsha Islands (Macclesfield Bank), Huangyan Island (Scarborough Shoal), and Nansha Islands (Spratly Islands). The total area of coral reefs in the SCS is approximately 37,935 km², accounting for 5% of the world’s coral reef area (Wang et al. 2014). Among the aforementioned five islands, Xisha Islands (15°46′N−17°08′N, 111°11′E−112°54′E) are located in the center of the SCS’s coral reef system as shown in Figure 1, and are the archipelago with the largest land area in the SCS.

Figure 1. (a) Coral reef distribution of the South China Sea (SCS, the location of Xisha Islands was derived based on results of this study, and the location of other reef areas was obtained from two resources introduced in section 2.2.2), and the location of Xisha Islands in the SCS is presented in dashed box. (b) Enlarged view of Xisha Islands, showing the spatial distribution of coral reefs based on the inventory established in this study.

2.2 Data resources

2.2.1 Satellite datasets

The Operational Land Imager (OLI) onboard the Landsat-8 satellite collects images with a 16-day repeat cycle and a 30-m resolution for seven spectral bands from visible to short-wave infrared (433 nm−2300 nm). Landsat-8 OLI captures data with improved radiometric precision over a 12-bit dynamic range (USGS 2019). Its improved signal to noise ratio (SNR) performance enables improved characterization of land/ocean state and condition, leading to better discrimination between coral reefs and water. Besides, it has been proved that Landsat-8 OLI can supply reliable remote sensing reflectance (Rrs) products, with improved geolocation accuracy (Wang et al. 2019), across shallow coral reefs and optically deep waters, allowing for geomorphic/benthic-habitat mapping and water column optical properties exploration (Hedley et al. 2018; Pahlevan et al. 2017; Wei et al. 2018).

The Xisha Islands are covered by four tiles of Landsat-8 OLI images based on the World Reference System-2 (WRS-2) path/row system (path/row: 122/48, 122/49, 123/48, 123/49). Focusing on long-term coral habitat change detection of Xisha Islands, a total of 763 Landsat-8 OLI images from 26 April 2013 to 31 December 2021 were used in this study, and the Level-1 scenes were downloaded from Earth Explorer of the United States Geological Survey (USGS) (https://earthexplorer.usgs.gov/).

2.2.2 Global datasets of coral reefs

In this study, the extent of coral reefs as well as lagoon and land area/location in the Xisha Islands were explored. For comparison purpose, the global inventory of coral reefs published by the United Nations Environment Programme World Conservation Monitoring Centre (UNEP-WCMC: https://data.unep-wcmc.org/datasets/1) (UNEP-WCMC, WorldFish Centre, WRI, and TNC 2021) was obtained. It is the most comprehensive global dataset of distribution of coral reefs in tropical and subtropical areas to date, serving as a baseline map for more detailed work in the future. Yet, this dataset was established based on Landsat series data before 2010, and reef distribution of Xisha Islands is not validated (VERIF flag: Not Reported). Thus, the administration boundaries of cities and regions in China provided by Resource and Environment Science and Data Centre (RESDC) were also downloaded (https://www.resdc.cn/) to provide more reference information. This is an official dataset containing border information of the administrative divisions of all provinces and cities in China and could be regarded as reliable geographic information data.

Moreover, as a public coral reef habitat mapping reference, the Allen Coral Atlas provides detailed classification results as an open resource (https://www.allencoralatlas.org/atlas/) for habitat map creation, dynamic monitoring, and is convenient for seeking regionally comparable data. The Allen Coral Atlas primarily relies on high-resolution satellite (e.g. WorldView-2, WorldView-3) imagery, which can capture detailed features of coral reefs, like coral health, reef structure, etc. Satellite data are made available after extensive processing with Global Discovery and Conservation Science (GDCS) algorithms and further transferred into habitat maps with 5 m resolution based on machine learning techniques and object-based analysis (Kennedy et al. 2021; Lee, Carder, and Arnone 2002; Lee, Weidemann, and Arnone 2013; Lee et al. 1999; Li et al. 2019; Lyons et al. 2020; Roelfsema et al. 2018, 2021). Quality assurance of this Atlas typically relies on ground truth data collection, expert review, and cross-validation with other data sources. Though not all results have been verified yet, more actual feedbacks by Atlas users/experts and repeated deliberation and misclassifications identification by Atlas team members will consistently improve the accuracy of output products. Presently, the overall accuracy of habitat mapping provided by Allen Coral Atlas at Qilianyu Islands is 80.6%, among which the recognition accuracy of coral reefs has reached to 88% (Tian, Zhu, and Han 2020).

2.2.3 Field survey data

A field survey was conducted on the Qilianyu Islands on 22 May 2023 and 24 May 2023. Water depth (including water quality) was obtained during fieldwork using the Multi-parameter Water Quality Monitor-OT7. A total of 27 survey sites were collected, each of which had the Global Positioning System (GPS) coordinates and water depth. In light of inevitable influence of tide in the real-time water depth measurement, tide correction has been implemented after data acquisition. The tide values were extracted and analyzed from the nearest tide gauge station located in Yongxing Island. Both the survey point locations and tide gauge station are shown in Figure 2.

Figure 2. Location map of the field survey points and tide gauge station. Field investigation conducted on 22 May 2023 and 24 May 2023 is depicted with yellow hollow triangle and orange solid square. Tide gauge station location is expressed with red pentagon in thumbnail. The background satellite image shown is a Landsat-8 OLI image acquired on 17 May 2023.

2.2.4 Environmental variable products

Monthly global 5 km satellite coral bleaching heat stress degree heating week (DHW) products are provided by NOAA Satellite and Information Service (https://coralreefwatch.noaa.gov/product/index.php). DHW generally range between 0 and 20°C-weeks, which can not only reflect accumulated heat stress over the most recent 12 weeks, but also indicate possible coral bleaching or death (McClanahan et al. 2020). DHW is closely related to the intensity and prolonged periods of coral bleaching. Coral bleaching occurs when DHW value reaches 4°C-weeks. By the time DHW value hits 8°C-weeks, severe, widespread bleaching is probable and significant mortality can be expected (Glynn and D’Croz 1990; Skirving et al. 2020).

Besides, taking the limited coverage of in-situ data into account, the Navigation Charts (v2020.6) of the Xisha Islands provided by China Navy Hydrographic Office was collected as well. The Navigation Charts was compiled based on historical voyage measured data and universal hydrographic data model (S-100), aiming at providing essential information such as water depth and terrain during vessel navigation. Therefore, it was regarded as a reference of water depth for the verification of the water depth calculation scheme in the following sections.

2.3 Development of mapping framework

The mapping framework applied in this study could be summarized in four parts: data pre-processing, probabilistic inundation mapping, geomorphic zonation/benthic habitat mapping, and data validation. A detailed schematic of each experimental procedure has been connected and presented in Figure 3. More specific content of each step was introduced below.

Figure 3. Schematic of all steps used to develop the mapping framework.

2.3.1 Data pre-processing

Data pre-processing procedures of Landsat-8 OLI images were divided into three sub-steps: (1) Atmospheric correction of Landsat-8 OLI images to retrieve the spectral remote-sensing reflectance (Rrs; sr⁻¹) from the original digital number (DN) was implemented in SeaWiFS Data Analysis System (SeaDAS, v7.5.3), which is an open-source software package distributed by National Aeronautics and Space Administration (NASA) and is widely applicable in processing, analyzing, quality control, and displaying of ocean color measured data from a host of satellite-based multi-spectral radiometers. (2) For removing the influence of irregular cloud, cloud pollution area was yield based on Quality Assessment (QA) bands provided by Landsat-8 OLI data. The final determined cloud pollution area was dilated with a radius of 3 pixels centered on each pixels identified as “cloudy” to prevent potential omission and false detection in the QA bands. (3) Considering each segmentation of the Xisha Islands has separate reef characteristics and ambient water environment, the dataset of the entire Xisha Islands was split into isolated image chips in order to reduce disturbance and save the subsequent computing time as well.

2.3.2 Probabilistic inundation mapping and reef extent extraction

Probabilistic inundation maps for coastal terrestrial ecosystems are generally used to identify flood prone areas (Alfonso, Mukolwe, and Di Baldassarre 2016; Maranzoni, D’Oria, and Mazzoleni 2022; Sarhadi, Soltani, and Modarres 2012). In a recent study, Dong et al. (2019) applied inundation mapping method based on time-series Landsat-8 OLI images and derived inundation frequency of coral reefs in the Nansha Islands (Spratly Islands), which could serve as an indirect depiction of bathymetry and further generate geomorphic zonation.

Besides, assuming that geomorphology of coral reefs is relatively steady within a relative short period when apart from human influences, the distinctions of reflectance between coral reefs and surrounding water mainly ascribes to water depth and sensor bands: (1) Admittedly, the varied water depth of coral reefs area is generally affected by occasional and transient phenomenon (e.g. typhoon and gyres) and periodic occurring tidal action. However, from the perspective of spatiotemporal average of extensive remote sensing images, the ephemeral constitutes can be receded, while the regular constitutes will be consolidated by the time-series manipulation. (2) Referring to influence of sensor bands, it is common sense that different bands have different penetrating depth, which is fundamental in distinguishing coral reef and its surrounding water and detecting geomorphic features. Since the light attenuation degree in water is specific to band, the longer the band, the shallower the penetration depth. As for Landsat-8 OLI, containing four visible bands (Band1: 430–450 nm, Band2: 450–515 nm, Band3: 525–600 nm, Band4: 630–680 nm), one near-infrared (NIR) band (Band5: 845–885 nm), and two shortwave-infrared (SWIR) bands (Band6: 1560–1660 nm, Band7: 2100–2300 nm), the penetrating depth varies from ~30 m (Band2 in clear water) to 0 m (SWIR bands are unable to go through water surface in the presence of strong absorption) (Jupp 1988). Accordingly, the area of inundation decreases with increasing bands and the depth range considered in the present work is up to 30 m.

Based on the aforementioned facts, probabilistic inundation mapping should be implemented at time-series image chips containing isolated reef (results in data pre-processing step) to mitigate the influence of optical properties particular to each segmentation. Given the position-invariant feature and adequate accuracy of cloud detection, accumulation of temporal sequence results of probabilistic inundation (PI) can not only delineate the inundation extent (boundary of remotely measurable area) of coral reefs, but also reflect their geomorphic characteristic. Higher the PI value, more frequent perceived by optical remote sensing sensors and shallower water depth at corresponding pixels. Inversely, pixels with lower PI value are presumed to be deeper part. Concrete PI value calculation can be described as the following steps:

1. Classification with Fuzzy ISODATA (FISODATA) method. FISODATA is an iterative clustering algorithm which has a capability of self-organizing by splitting and merging clusters (Bezdek 1980). Compared to conventional hard clustering methods, FISODATA gives the flexibility of mapping a pixel to one or more clusters rather than assigning each pixel to one explicit class, and a membership or affinity score is meanwhile introduced to determine final attribution of each pixel (Sreevalsan-Nair 2020; Zhao et al. 2022). FISODATA method has been confirmed to be competent in segmentation of multispectral images (Lei et al. 2018; Yang et al. 2015). In addition, it requires rarely priori input which is likely to be the case for remote coral reef areas. Therefore, FISODATA is primarily performed here to distinguish above-water and underwater parts (relative to the penetration depth of the band) of coral reefs at each Landsat-8 OLI band, and generate binary classification results for next accumulation step. Notably, affected by waves, fog, and sun glint, inappropriate binary classification consequence is inevitable. Hence, all of the binary image chips were visually inspected to ensure that only the correct results were conserved for accumulation.

2. PI value calculation. The time-series binary images chips from OLI band 1-6 were first accumulated for each segmentation of the Xisha Islands respectively. Then the results of each band were accumulated in union (accumulation and converted to probability between 0 and 1). For the purpose of making the calculation results comparable among separate segmentations, as well as different pixels, the PI value calculation results at each pixel were acquired by dividing the number of non-cloud observations at corresponding location. More detailed information about PI value calculation could be found in the work by Dong et al. (2019).

As mentioned before, PI results could be considered as a substitute of water depth. In order to verify this statement, the average of PI values within 3 × 3-pixel window at each survey point were derived for comparison purpose in subsequent analysis. To supplement this viewpoint, 92 points were manually selected and their depth values in the Navigational Charts of the Xisha Islands were extracted to make a further comparison between PI value and channel depth. In light of the comparability of data, the average of PI values within a neighbor window of 3 × 3 pixels centered at each selected point was applied likewise.

Take a further step, in the light of position-invariant characteristic of reefs within a comparatively short period and satisfactory geometric accuracy of Landsat-8 OLI images (~18 m) (Liu et al. 2019), reef extent at time-average level can be identified with accumulation results of binary images generated from Rrs retrieved by Landsat-8 OLI Bands 1–6. After we acquired PI results, the extent of individual reef as well as the boundary of land or lagoon were extracted afterward with traditional mathematical morphology operation including threshold segmentation, edge detection algorithms (here we applied canny operator), and dilation and erosion algorithm. All these procedures were conducted with built-in functions in MATLAB (r2020a). An example (Zhongjian Island) of extent recognition procedure was interpreted in Fig. S2 in Supplementary Material.

2.3.3 Water column correction

Sea surface reflectance after atmospheric correction of each segmentation image part was calculated first, and PI value was subsequently derived and applied as a water depth substitution index. However, the interpretation of benthic information was still impeded by the fact that bottom reflectance signal was difficult to be distinguished from water depth fluctuations. Therefore, it is necessary to remove the influence of water depth variation before we realize accurate benthic classification.

Spectrum varies with water depth, leading to the universal phenomenon that same objects are with different spectrums or same spectrums stand for different objects. In order to better extract bottom features, the influence of water depth need to be removed first. Generally, the optical properties of water column can be considered to be vertical and horizontal homogeneous in certain clear shallow water area, while the depth of the water is highly variable. Based on this assumption, Lyzenga (1978, 1981) proposed a water column correction method to combine the information in different spectral bands to produce a depth-invariant index ($DI{I}_{ij}$) as an alternative of the bottom type. The specific depth-invariant index retrieval procedure can be summarized as the following steps:

1. Logarithm transformation. The method was developed on the basis of a simple water reflectance model, which assumed that light attenuation exhibited an exponential decrease with increasing depth. Naturally, the relationship between the transformed reflectance and water depth could be approximately linearized with logarithm transformation (Eq. 1, $L_{si}$ is the deep-water reflectance at band i, $L_i\hspace{0.17em}$is the shallow-water reflectance at band i).

   (1) ${\mathrm{X}}_{\mathrm{i}}=\mathrm{l}\mathrm{n}({\mathrm{L}}_{\mathrm{i}}\text{-}{\mathrm{L}}_{\mathrm{s}\mathrm{i}}),$(1)

2. Attenuation coefficient ratio. A proxy of the attenuation coefficient ratio between band i and band j ($K_{ij}$) was determined as regression slope of a bi-plot between transformed reflectance in the two bands ($X_i$ and $X_j$), and it is a constant value for any benthic type. In this study, sandy area (sand was always displayed as highlighted part in Landsat-8 OLI image which was easy to be identified) was selected as uniform substratum with variable water depth to calculate $K_{ij}$ (Eq. 2 and Eq. 3). In Eq.2 and Eq.3, $K_i$ or $K_j$ stands for irradiance attenuation coefficient of the water at band i/j, ${\sigma }_{ii}$ or ${\sigma }_{jj}$ represents reflectance variance of band i/j.

   (2) ${\mathrm{K}}_{\mathrm{i}\mathrm{j}}=\frac{{\mathrm{K}}_{\mathrm{i}}}{{\mathrm{K}}_{\mathrm{j}}}=\mathrm{a}+\sqrt{{\mathrm{a}}^2+1},$(2)
   (3) $\mathrm{a}=\frac{{\sigma }_{\mathrm{j}\mathrm{j}}-{\sigma }_{\mathrm{i}\mathrm{i}}}{2{\sigma }_{\mathrm{i}\mathrm{j}}}\mathrm{a}\mathrm{n}\mathrm{d}{\sigma }_{\mathrm{i}\mathrm{j}}=\overline{{\mathrm{X}}_{\mathrm{i}}{\mathrm{X}}_{\mathrm{j}}}-\overline{{\mathrm{X}}_{\mathrm{i}}}\overline{{\mathrm{X}}_{\mathrm{j}}},$(3)

3. Depth-invariant index. Finally, the depth-invariant index ($DI{I}_{ij}$), which is unique to two bands, can be generated according to Eq. 4.

   (4) $\mathrm{D}\mathrm{I}{\mathrm{I}}_{\mathrm{i}\mathrm{j}}={\mathrm{X}}_{\mathrm{i}}\text{-}{\mathrm{K}}_{\mathrm{i}\mathrm{j}}{\mathrm{X}}_{\mathrm{j}},$(4)

After water column correction, $DI{I}_{12}$, $DI{I}_{13}$, $DI{I}_{14}$, $DI{I}_{23}$, $DI{I}_{24}$, $DI{I}_{34}$ were produced and can be employed as input variables in the subsequent benthic habitat classification process combined with probabilistic inundation results.

2.3.4 Coral reef geomorphic zonation and benthic habitat classification

The establishment of an appropriate classification system is critical in transforming remote sensing observations into user-friendly spatial information. In this study, a three hierarchical classification scheme was first applied for reef extent detection followed by specific geomorphic zonation and benthic cover type classification: (1) Reef Extent; (2) Geomorphic Zonation; and (3) Benthic Cover Type. Hierarchical structure, mapping category, detail description, and mapping rules are summarized in Table 1. Reef extent was identified with PI product which reveal the water depth variation combined with OLI band 1–6 exploration range. Afterwards, incorporated with PI results as well as slope derived from PI results, FISODATA method introduced in Section 2.3.2 was applied to map geomorphic zones under a protocol provided by Roelfsema et al. 2013, 2018), which comprehensively considers the formation mechanism of coral reef, the morphology of coral reef, and the interpretability degree of remote sensing images. FISODATA requires a spatial dataset as input to generate segmented groups of pixels with resemble feature (here mainly refers to PI results and slope), followed by labeling using a membership rule shown in Table 1. Similarly, dominant benthic cover type was also classified with FISODATA method while based on sub-surface reflectance after water column correction described in Section 2.3.3 and PI results. The classification criterion of benthic cover type was originated from previous mentioned protocol with appropriate modification considering local knowledge and data characteristic (Jiang et al. 2018; Li 2021; Zuo et al. 2017), and finally type labels assigned included Coral/Algae, Seagrass, Microalgal Mats, Rock, Rubble, and Sand. Taking Bei Reef as an example, the geomorphic zonation mapping process and the benthic habitat mapping process were described in Fig. S3 and Fig. S4, separately.

Table 1. Hierarchical structure, description, and mapping rules (PI = probabilistic inundation value).

Hierarchical levels	Mapping category	Detail description	Mapping rules
Reef extent	Land	Exposed above water	PI = 1
	Reef Shallow	Emergent at low tide	0 < PI < 1
	Deep Water	Submerged even at low tide	PI = 0 in open area
Geomorphic zonation	Reef Slope	Steep sloping area extending seaward	Based on classification results and adjacent to Deep Water
	Reef Front	Upper part of the seaward slope	Based on classification results and adjacent to Reef Slope
	Reef Crest	The boundary between the Reef Flat and Reef Front	Based on classification results and adjacent to Reef Front
	Outer Reef Flat	Adjacent to the Reef Crest and is a shallow carbonate platform	Based on classification results and adjacent to Reef Crest
	Inner Reef Flat	Gently sloping area adjacent to the Outer Reef Crest	Based on classification results and adjacent to Outer Reef Flat and Lagoon
	Back Reef Slope	Interior area behind the Inner Reef Flat and is sediment-dominated	Based on classification results and between Inner Reef Flat and Lagoon
	Shallow Lagoon	Semi to fully-enclosed and sediment-dominated lagoon area with water depth shallower than 5 m approx.	Based on classification results and adjacent to Back Reef Slope
	Lagoon	Semi to fully-enclose by the Reef Flat with variable water depth deeper than 5 m approx.	Based on classification results and adjacent to Shallow Lagoon
	Patch Reef	Isolated and small reef outcrop within lagoon	Based on classification results and distinguishable in Lagoon area
Benthic cover type	Coral/Algae	Commonly associated with Reef Front, Reef Slope, Patch Reef and Outer Reef Flat	Based on classification results and Remote Sensing Images
	Seagrass	Commonly associated with Shallow Lagoon and Back Reef Slope	Based on classification results and Remote Sensing Images
	Microalgal Mats	Commonly associated with sandy area	Based on classification results and Remote Sensing Images
	Rock	Commonly associated with Reef Crest and Outer Reef Flat with pavement appearance	Based on classification results and Remote Sensing Images
	Rubble	Commonly associated with Reef Crest and Reef Front	Based on classification results and Remote Sensing Images
	Sand	Commonly associated with Inner Reef Flat, Back Reef Slope and Lagoon	Based on classification results and Remote Sensing Images

Since different islands have diverse properties, rules varied in different islands in threshold value for a dominant benthic cover type. Generally speaking, a sensible classification system should make a sound ecological sense, like spectrally separable between different benthic types (Xu and Zhao 2014). However, restricted by spatial resolution and nonnegligible noise in remote sensing data, the results of geomorphic zonation and benthic type classification were not precise, which should be regarded as an approximation of reality for reference purpose.

It is also worth noting that the geomorphic zonation and benthic habitat data provided by Allen Coral Atlas of the Xisha Islands as validation reference were retrieved in 2020; given the consideration that the benthic habitat types were relatively susceptible to environmental changes or artificial interference, the geomorphic zonation and benthic habitat results applied in comparison were also only built with Landsat-8 OLI images in 2020. Please refer to Table S3 for more detailed information about Landsat-8 OLI images applied to retrieve geomorphic zonation and benthic habitat results in 2020. Moreover, in contemplation of maintaining consistent spatial resolution with OLI images, the geomorphic and benthic data downloaded from Allen Coral Atlas in polygon format with geometry unit information were resampled and remapped as well.

3 Results

3.1 Probabilistic inundation mapping of the Xisha Islands

Probabilistic inundation maps of coral reefs in the Xisha Islands were shown in Figure 4. At a global scale, the PI of entire coral reefs from low to high values was shown with a color gradient from dark blue to light blue, reflecting the actual shifting in bathymetry to some extent. To display more detailed inundation situation, individual scale-up exploration was performed with each reef and was shown with a distinct designed zoom-in gradient color scheme (mainly from light blue to deep orange with four transition color sections). Permanently exposed areas (PI = 1), such as land and islet, which can always be detected in Bands 1–6, were shown with deep orange. The PI value decreased gradually from these permanently exposed area to deep water, and their spatial variation patterns were similar to the underwater topography.

Figure 4. Probabilistic inundation map of the Xisha Islands (PI stands for probabilistic inundation value).

In order to assess the consistency between PI value and water depth, in-situ data of Qilianyu Islands collected on 22/24 May 2023 were adopted to compare with PI values centered at each sample point. It is of worth attention that according to characteristics of PI calculation procedure, PI results indicated the spatial and temporal averaged immersion situation, which smooth the ephemeral constitutes and consolidated the periodic elements. Thus, probabilistic inundation map, to some extent, reflects the mean sea level in the Xisha Islands. Simultaneously, in-situ water depth data have been tidal-corrected and could be considered as mean sea level. Therefore, PI value and in-situ water depth data were comparable and the comparison results were presented in Figure 5(a). Moreover, taking into account that spatial coverage of measured data was limited, the correlation between PI value and channel depth from Navigational Charts were further analyzed and depicted in Figure 5(b) (for more details, please refer to Fig. S1 in Supplementary Material). As can be seen from Figure 5(a), an obvious linear regression trend was determined between PI value and in-situ data with the R² of 0.91 and the root mean square error (RMSE) of 1.36 m. Correspondingly, in Figure 5(b), a linear relationship was also discovered between PI and channel depth with the R² of 0.65 and RMSE of 4.16 m, though the fitting points are not densely concentrated. This could be partly attributed to the accuracy of channel depth and the deviation introduced by manually selecting verification points, as well as affected by the distribution of selected locations since the relatively dispersed points might suppress regression effects. Overall, both the correlation between PI and in-situ water depth around Qilianyu Islands and that between PI and channel depth in Xisha Islands showed satisfied linear relation in good consistency. The above comparison results also exhibited resemblance to the results by Dong et al. (2019), which further demonstrated the transferability of PI mapping method and indicated that the PI derived from time-series Landsat-8 OLI images dependably delineated the variation of water depth in the Xisha Islands.

Figure 5. (a) Scatterplot and fitted regression line for the relationship between PI value and in-situ water depth (27 pairs); (b) Scatterplot and fitted regression line (black solid line) for the relationship between PI value and channel depth extracted from the Navigational Charts (92 pairs). Three dashed lines (d1–3) were the results of Dong et al. (2019) and were shown here for the purpose of comparison.

3.2 Reef extent mapping

Reef extent, including boundary of land or lagoon, was extracted on the basis of PI results (Figure 6). According to reef extent mapping, an inventory containing the abundance and detected coral reef names both in Chinese and English was also given in Table S1 in the Supplementary Material. The geographical names of each detected coral reef wer referred to recorded historical information (Geographical Names Committee of Guangdong Province 1987). In the sake of convenience, Chinese names will be employed in the henceforth description. Notably, careful visual inspection was applied to inspect the location and existence of coral reefs in the Xisha Islands, and the inventory only retained reefs that are located in water depth shallower than 30 m within the penetration depth of optical remote sensing. A total of 38 isolated islands, sands, banks, reefs, and shoals were confirmed from the time-series Landsat-8 OLI images of the Xisha Islands during the period 2013–2021. According to the results, coral reefs in the Xisha Islands cover an area (within the detection range of Landsat-8 OLI Band1–6) about 560 km². The distribution of coral reefs in the Xisha Islands can be divided into two parts: Yongle Islands (eastern part) and Xuande Islands (western part), extending from Bei Reef (~17.13°N) to Songtao Bank (~15.72°N), and from Zhongjian Island (~111.18°E) to Xidu Bank (~112.9°E) (Figure 1b).

Figure 6. Comparison of extracted boundaries of coral reefs (including lagoon and land above water) in this research and the RESDC/UNEP-WCMC dataset. The RGB images presented the boundaries validation results of Bei Reef (red box), Qilianyu Islands (orange box), Yuzhuo Reef (blue box), and Zhongjian Island (green box) at two different years, respectively.

Additionally, collected reef extent products (land/lagoon area was also incorporated) were further compared with another two datasets of coral reefs: the UNEP-WCMC dataset and the RESDC dataset (Figure 6). Generally speaking, the reef extent exhibited resemble distribution spatially, when several existed obvious deviations, like the extent of Dongdao Reefs and Panshiyu Island between extraction results of this study and of the RESDC dataset. The deviation of results can be primarily attributed to the differences in data sources and data acquisition time. Additionally, derived reef extent in this research displayed more specified geomorphic information comprising lagoon, land, and hidden reefs without condition in which the outlines were discontinuous or jagged. In order to reflect intuitive and reliable comparison, we chose four reef sites (Bei Reef, Qilianyu Islands, Yuzhuo Reef, and Zhongjian Island) and corresponding cloud-free RGB images at two different years (Bei Reef: 2013/2021, Qilianyu Islands: 2014/2021, Yuzhuo Reef: 2013/2021, Zhongjian Island: 2013/2021). RGB images combined with reef boundaries were also exhibited in Figure 6, which further demonstrated that reef boundary extraction results of this study were more in line with the actual boundary conditions whether in an earlier year or a more recent year. Based on the above facts, the reef extent mapping provided in this study can be regarded as a practical implement in coral reef distribution status identification.

3.3 Geomorphic zonation mapping and area statistics of the Xisha Islands

Since PI could effectively reflect the variation of water depth, which also implied the underwater topography, PI and the slope calculated from PI were applied as input of FISODATA to produce geomorphic zonation map. Considering the unavoidable classification errors (Phinn, Roelfsema, and Mumby 2012; Roelfsema et al. 2018), the ultimate results presented have been visually corrected to ensure the accuracy. Finally, the geomorphic zonation was composed of 10 components: Reef Slope, Reef Front, Reef Crest, Outer Reef Flat, Inner Reef Flat, Back Reef Slope, Lagoon, Shallow Lagoon, Patch Reef, and Land (Figure 7). Besides, the concrete attributes of each geomorphic zone category were summarized in Table 2.

Figure 7. Geomorphic zonation map of the Xisha Islands (PI stands for probabilistic inundation value and the color scheme refers to the Allen coral Atlas).

Table 2. Area of major coral reefs of the Xisha Islands (km^2.).

Name	Exposed Area	Reef Flat Area	Sub-tidal Reef Area	Lagoon Area
Bei Reef	-	27.73	1.88	10.63
Jinyin Island	0.38	8.91	3.94	-
Ganquan Island	0.36	0.20	0.94	0.35
Shanhu Island	0.31	3.03	1.88	2.00
Quanfu Island	0.03	2.80	4.25	4.60
Yin Island	0.01	6.99	3.41	2.70
Jinqing Island	0.22	16.55	5.62	36.95
Chenhang Island	0.72	1.21	1.15	1.16
Lingyang Reef	0.02	8.83	2.26	5.58
Yuzhuo Reef	-	30.85	8.98	8.74
Huaguang Reef	-	50.46	65.90	111.91
Panshiyu Island	0.06	16.86	8.37	6.90
Zhongjian Island	1.28	4.25	1.25	1.20
Qilianyu Reefs	1.23	23.18	31.82	16.10
Yongxing Island	3.46	3.07	21.16	20.97
Dong Island	1.98	1.92	48.16	8.62
Langhua Reef	-	37.23	20.11	27.88

Exposed Area represents land area; Reef Flat Area represents Inner Reef Flat area and Outer Reef Flat area; Sub-tidal Reef Area represents Reef Slope area, Reef Front area and Reef Crest area; Lagoon Area represents Back Reef Slope area, Lagoon area, Shallow Lagoon area, and Patch Reef area.

Quantitative comparison results between geomorphic zonation classification in this study and that of Allen Coral Atlas were summarized in Table S2. The overall accuracy varied between 49% (Zhongjian Island) and 81% (Panshiyu Island) with 72% average accuracy of the entire Xisha Islands (Table S2). The relative low classification accuracy of Zhongjian Island was attribute to images of this area affected by clouds and their shadows frequently and seriously. It led to discontinuous and contrast-lacking of PI results at this area, introducing to ultimate misclassification. Confronted with reefs with larger reef area, reefs with small area were more susceptible to changing environment, resulting in poor classification results. The overall geomorphic zonation classification accuracy not reaching 70% was concentrated in reefs with smaller area (Dong Island, Zhongjian Island, Shanhu Island, Quanfu Island, Jinqing Island, Chenhang Island (for detailed information, please refer to Table S2)). According to the geomorphic zonation results, exposed area (land), reef flat area (outer reef flat, inner reef flat), sub-tidal reef area (reef slope, reef front, reef crest), and lagoon area (back reef slope, lagoon area, shallow lagoon, and patch reef) were calculated and summarized in Table 2. It can be seen from Table 2 that Yongle Atoll, as a relative complete reef system, owned the largest reef area close to 130 km² (only refers to the part that can be detected by optical remote sensing) and contained the most comprehensive geomorphology. Besides, Huaguang Reef was an independent island with the largest lagoon area (~115 km²), while Yongxing Island was an independent island with the largest land area (~3.5 km²).

3.4 Benthic habitat mapping

Similar to geomorphic zonation mapping procedure mentioned above, here we combined PI, DII, and four Rrs (443, 482, 561, and 655 nm) as input of FISODATA to generate benthic habitat map. Eventually, the benthic habitat results derived from FISODATA were regrouped and assigned with eight labels: Land, Lagoon/Deep Water, Coral, Seagrass, Microalgal Mats, Rock, Rubble, and Sand (Figure 8). In contrast to geomorphic zonation mapping, the essential difference could be seen from Figures 7 and 8 that benthic habitat classification results did not cover whole area under the influence of cloud blocking or data noise. For most reefs, several (often 3–5) images of similar weather conditions could be found, which was convenient to integrate to make up for local data lacking. However, there were still some reefs that did not have clear images or the weather condition between different images were of comparatively great difference interfering data fusion, leading to the absence of data shown in Figure 8.

Figure 8. Benthic habitat map of the Xisha Islands (PI stands for probabilistic inundation value and the color scheme refers to the Allen coral Atlas).

In order to make a comparison with benthic habitat data provided by Allen Coral Atlas of the Xisha Islands, the benthic habitat results presented in Figure 8 were also only built with Landsat-8 OLI images in 2020. Preciseness of the benthic habitat classification results compared to that of Allen Coral Atlas varied between 76% (Bei Island) and 94% (Yinyu Island/Yinyu Islet) with 86% average overall accuracy of the entire Xisha Islands (Table S4). It was also important to note that only using PI with Rrs as input of FISODATA, without consideration of water depth variation, may not get benthic classification results with acceptable accuracy. To illustrate the necessity of introducing DII during benthic classification procedure, we compared two classification results with or without DII to the Allen Coral Atlas at Bei Reef, and the confusion matrix was shown in Tables 3 and 4. Clearly, it improved classification accuracy when the effects of variation in water depth were taken into account.

Table 3. Confusion matrix for benthic cover type classification results at Bei Reef without taking water depth variation into consideration (I: Land; II: Coral/Algae; III: Seagrass; IV: microalgal Mats; V: Rock; VI: Rubble; VII: sand).

		This research
Bei Reef		I	II	III	IV	V	VI	VII
Allen coral atlas data	I	-	-	-	-	-	-	-
	II	-	3714	75	131	1897	13	43
	III	-	136	1023	246	136	71	178
	IV	-	363	19	1987	43	103	31
	V	-	1141	3	11	7849	1784	98
	VI	-	451	19	1063	2131	9787	526
	VII	-	25	0	56	332	584	2673
	False Negative	-	0.36	0.01	0.43	0.37	0.21	0.25
	True Positive	-	0.64	0.90	0.57	0.63	0.79	0.75
	Accuracy	0.70

Table 4. Confusion matrix for benthic cover type classification results at Bei Reef when taking water depth variation into consideration (I: Land; II: Coral/Algae; III: Seagrass; IV: microalgal Mats; V: Rock; VI: Rubble; VII: sand).

		This research
Bei Reef		I	II	III	IV	V	VI	VII
Allen coral atlas data	I	-	-	-	-	-	-	-
	II	-	4819	0	0	1871	13	29
	III	-	21	1139	0	141	12	225
	IV	-	30	0	2504	4	124	0
	V	-	833	0	0	8041	1566	98
	VI	-	103	0	914	2244	10030	420
	VII	-	24	0	76	87	597	2777
	False Negative	-	0.17	0	0.28	0.35	0.19	0.22
	True Positive	-	0.83	1	0.72	0.65	0.81	0.78
	Accuracy	0.76

In addition, the density distribution of PI value corresponding to each benthic habitat type at each individual reef as well as the entire Xisha Islands was further explored (Figure 9). As for reefs containing land area (Xuande Reefs, Dong Island, Yongle Atoll, Panshiyu Island, and Zhongjian Island), the density distribution of PI value about land area mainly concentrated in the part with PI value greater than 0.8, indicating that the long-term average exposure degree (could be detected with optical remote sensing) of this region was relatively higher than other parts. This was consilient with objective conditions and further verified that PI value exploited as an indicator of water depth and employed to geomorphic zonation/benthic habitat classification was applicable. Another significant feature shown in Figure 9 was that despite coral distributed in various water depth, they were primarily concentrated in relative deep-water areas (PI < 0.2, the corresponding water depth was deeper than 15 m) with more fitting temperature, relatively stable environment (not susceptible to waves or extreme weather), and sufficient settling organic matter. As for other benthic habitat types, they were evenly spread in different water depth without apparent relationship with water depth.

Figure 9. Probability density distribution map of relationship between benthic habitat type and probabilistic inundation (PI) of the Xisha Islands.

4 Discussion

4.1 Application of probabilistic inundation mapping and its potential limitations

In consideration of hard in-situ investigation of the Xisha Islands, deficiency or even absence of field data in the spatio-temporal domain was a common situation. In this study, we employed a simple method using time series Landsat-8 images from 2013 to 2021 to build up probabilistic inundation map within penetration depth of optical remote sensing. Correlation analysis was implemented between PI results and in-situ water depth data collected at Qilianyu Islands, as well as that between PI results and channel depth derived from Navigational Charts in the Xisha Islands, to validate and determine whether the variation of PI was consistent with water depth fluctuation of coral reefs. The comparison results turned out that PI derived from time-series Landsat-8 OLI images authentically depicted the variation of water depth around the coral reef in the Xisha Islands.

The method of PI mapping of coral reefs introduced hereinbefore was on the premise of several hypothesis: (1) high-quality Landsat-8 OLI images with sufficient geometric accuracy were available, allowing for the direct accumulation of the extraction results of PI from multiple images; (2) besides, sufficient number of Landsat-8 OLI images obtained at different water environment (various water levels affected by weather or tides) were possible to be collected, enabling the production of quantitative and continuous time series cumulative result which can represent the geomorphology of individual coral reef and also reflect the full view of the Xisha Islands; (3) the topography was relatively stable during the observation period since the transformation of coral reefs (growth or degradation) is a long-term process; (4) the extent area, which was extracted from PI results delineating boundaries between coral reefs and deep water, could exactly represent the inundation under objective conditions within the penetration depth of optical wavelengths.

It is noteworthy that the exactness of PI may be affected by the stability of sands. According to Li et al. (2022), the geomorphology of some sands was not stable yet. Due to different positions in reef flat and types of sediments, each sand is subject to different hydrodynamic conditions; thus, the evolution process and stable state were inconsistent with each other (Perry et al. 2011). Therefore, despite the fact that PI derived in this study could reflect the geomorphology of the Xisha Islands, it was a spatiotemporal average result closely related to the time span and location of Landsat data used in the study. For regions where the geomorphology has not maturely developed at this stage, a longer period of monitoring and investigation of the surrounding natural environment are needed in the future. Admittedly, since data noise and inversion error of Rrs (influenced a lot by colored dissolved organic matter or total suspended matter) were unavoidable and uncertainties in data processing or misidentification of classification procedures would potentially introduce errors into the final results (Kennedy et al. 2021), PI accumulated with binary results of Rrs at each band using FISODATA method showed varied accuracy at different regions, while here we only considered the overall accuracy and regarded it as a reliable alternative geomorphology. More detailed study referring remote sensing data processing and classification method comparison should be explored in the future.

4.2 Qualitative and quantitative analysis of geomorphologic features of typical reefs

The geomorphic zonation depicted in Figure 7 was built with Landsat-8 OLI images from 2013 to 2021 since PI was long-term average value, while the validation geomorphic data provided by Allen Coral Atlas of the Xisha Islands were retrieved in 2020. Since the geomorphology formation of coral reefs is a long-term historical process, fundamental transformation will not occur in a short time without extreme external interference. Therefore, the geomorphic zonation of the Xisha Islands in 2020 explored in this research could be regarded as a stable characteristic during the whole study period. Coral reefs that have undergone destructive development or land reclamation should be considered separately, and these external factors were not included in this research since no obvious structural change was found during visual inspection.

Generally speaking, the classification results were predominantly consistent in the distribution at different geographical units, and it represented a common geomorphic feature that geomorphy of coral reefs of the Xisha Islands distributed in the NE-SW direction appear oblong or spindle-shaped, which was basically affected by the northeast and southwest monsoon in winter and summer, respectively (Xu et al. 2016). Hence most reefs have broad rims at the northeast and the southwest areas. Meantime, the rims at northeast part were comparatively broader than southwest part, as the impact of northeast monsoon has a stronger intensity and a longer duration than that of the southwest monsoon.

Among all the reefs in the Xisha Islands, there are seven typical individual reefs containing lagoon (Bei Reef, Qilianyu Reefs, Langhua Reef, Lingyang Reef, Yuzhuo Reef, Panshiyu Island, and Huaguang Reef). A qualitative and quantitative analysis was performed to explore the geomorphic characteristics of different reefs and to verify whether PI and derived slope could embody plausible geomorphic features. In contemplation of making PI results more intuitive, three-dimension structure and vertical view of reefs were depicted based on PI. Besides, profiles were selected in all seven reefs and PI on each profile was extracted with nearest neighbor interpolation to observe the variation of PI related with geomorphology (Figure 10). Enclosed (lagoon was inside of the reef flats) or semi-enclosed (lagoon was inside of reef flats with one or more entrances to deep water) reefs were easily to be extracted from PI results since mean water depth of lagoon was obviously deeper than that of reef flat part; hence, PI of lagoon should generally be larger than PI of reef flat. The distinction of PI between lagoons and reef flats was well reflected in Figure 10 a3, c3–g3, and PI of lagoons varied between 0 and 0.4 overall. Besides, in Figure 10 b3, PI of profile selected in Qilianyu Reefs also presented with satisfied discrepancy between land and reef flat. Land within the reef flat was exposed to air even at high water level; thus, PI of land area should exhibit high value, and PI = 1 at Xi Sand was in accordance with such fact. More detailed information can be derived from Figure 10 a1–g1 that PI could not only catch geomorphology characteristics with significant differences (land or lagoon), but also describe the topographic fluctuation within lagoon (patch reef) or the reef flat (morphological changes caused by sedimentation or waves).

Figure 10. Maps of three-dimension structure (a1–g1), vertical view (a2–g2), and selected profiles across the reefs (a3–g3) of seven reefs in the Xisha Islands based on probabilistic inundation (PI).

For the purpose of investigating the natural evolution of reefs, a quantitative analysis containing calculation of the total area (S), perimeter (P), area of lagoon (Sl), and area of reef flats (Sf) was accomplished from the PI results (Table 5). Moreover, the openness (O) and compactness (C) used to characterize the morphology of reefs were also determined according to the following definition (Liu 2001):

Table 5. Geomorphologic parameters of seven reefs in the Xisha Islands.

	P (km)	S (km^2)	Sl (km^2)	Sf (km^2)	O	C
Bei Reef	28.27	44.40	10.63	33.77	0.31	4.44
Qilianyu Reefs	42.00	34.01	1.61	32.40	0.05	3.19
Langhua Reef	39.69	70.55	27.88	42.67	0.65	4.73
Lingyang Reef	16.40	16.78	5.58	11.20	0.50	3.59
Yuzhuo Reef	33.61	43.23	8.74	34.50	0.25	4.02
Panshiyu Island	20.22	25.88	6.94	18.94	0.37	4.01
Huaguang Reef	64.01	182.69	111.91	70.78	1.58	5.99

P is the perimeter; S is the total area of reef; Sl is lagoon area; Sf is reef flat area; O is the openness; and C is the compactness.

(5) $O=\frac{Sl}{Sf}\hspace{0.17em},\hspace{0.17em}C=2\sqrt{\pi S/P}$(5)

The area derived here was simply pixel-based without consideration of distortion of remote sensing images; thus, it should be cognized as an approximation from an overhead view instead of real area taking the slope of terrain into account. As can be seen from Table 5 and Figure 10, the PI of lagoons was discovered to vary between 0.1 and 0.4 with lower value of compactness, while PI was normally lower than 0.1 with higher value of compactness. This phenomenon suggested that compactness may have a probable correlation with PI or water depth. As a comparison between lagoon area and reef flat area, openness assessed the development degree of reefs. Taking Huaguang Reef as an example; it has four entrances connecting lagoon with outer deep water, which accelerate the rate of water exchange, leading to be more vulnerable to the impact of environment variation (Tamborski et al. 2019). Despite owning strong connection with surrounding water and under richer nutritional conditions, Huaguang Reef is one of less developed reefs in the Xisha Islands with relatively narrow width of reef flats. However, on the contrary, patch reefs were prone to form in open lagoons as water fluxes might hasten the accumulation of coral secretions at lagoon bottom (Dong et al. 2019). Therefore, there exist plenty of patch reefs in relative open lagoons (e.g. Huaguang Reef, Langhua Reef), while patch reefs in highly developed reefs (e.g. Bei Reef) may gradually develop into a more nature state without external disturbance.

The reef flat area, lagoon area, land area, and total area of different reefs acquired in this study were further confirmed with similar previous works (Shen 2013; Xu 2012; Xu et al. 2016) and the comparison results were concluded in Tables 6 and 7. In Table 6, total area of Yongle Atoll was not documented because only part of lagoon area was explored in this study and most area of lagoon connected to deep water within Yongle Atoll was not considered, while it was all calculated in previous study. For other isolated reef containing lagoon, the enclosed lagoon areas of this study were consistent with that of Xu et al. (2016). Total area of Bei Reef, Yuzhuo Reef, Panshiyu Island, and Huaguang Reef derived in this study were 40.24, 48.57, 32.19, and 228.27 km², respectively, and that obtained by Xu et al. (2016) were 47.72, 48.96, 27.81, and 190.68 km², respectively. The differences between two research results were −7.48, −0.39, 4.38, and 37.59 km² and the mean deviation was 2.67%. The total area of each reef in this study was slightly different with previous consequence in view of different Landsat-8 OLI images applied introducing deviation substantially caused by tides or weather conditions. With respect to land area, the comparison between several research was accomplished and shown in Table 7. Evidently, the statistical land area incorporated in individual reefs was basically identical although researches were conducted in different years and in different measurements. The land area of Chenhang Island and Guangjin Island in this work was larger than preceding works (0.51, 0.38, and 0.18 km² larger than that by Xu (2012), Shen (2013), and Xu et al. (2016)) as the embankment bridging these two islands located on the same reef was also ascertained and measured in this study. In consonance with anterior studies, the geomorphic classification results of this paper based on time-series Landsat-8 OLI images effectively and precisely reflected the objective geomorphic information, especially the classification results of land and enclosed-lagoon showed marked differences with surrounding environment. The similarity of results in different researches also proved that the time span involved in this study was relatively short, during which the geomorphology of coral reefs was relatively steady and could be regarded as constant or with inconspicuous and undetectable variation. Therefore, the premise assumption for obtaining PI was valid. In contrast to previous studies, the area estimation in this study was based not on certain remote sensing image, but from time-series remote sensing images, implying that results obtained reflected the overall average state of all reefs, which were objective and referential without being limited by time and space.

Table 6. Comparison results of reef flat area, lagoon area, and total area between this study and previous study (Xu et al. 2016) (unit: km^2.).

	Reef flat area		Lagoon area		Total area
	(Xu et al. 2016)	This study	(Xu et al. 2016)	This study	(Xu et al. 2016)	This study
Bei Reef	30.92	27.73	10.96	10.63	47.72	40.24
Yongle Atoll	57.12	48.52	5.6	5.58	-	-
Yuzhuo Reef	34.06	30.85	8.75	8.74	48.96	48.57
Panshiyu Island	17.83	16.86	7.19	6.90	27.81	32.19
Huaguang Reef	63.07	50.46	113.62	111.91	190.68	228.27

Table 7. Comparison results of land area in different studies (unit: km^2.).

	(Xu 2012)	(Shen 2013)	(Xu et al. 2016)	This study
Jinyin Island	0.33	0.36	0.41	0.38
Ganquan Island	0.30	0.31	0.37	0.36
Shanhu Island	-	0.30	0.35	0.31
Quanfu Island	-	0.02	0.03	0.03
Yin Island	-	0.01	0.03	0.01
Jinqing Island	-	0.21	0.25	0.22
Chenhang Island/Guangjin Island	0.21	0.34	0.54	0.72
Lingyang Reef	0.42	0.01	0.01	0.02
Panshiyu Island	-	0.4	0.06	0.06

4.3 Trends of coral area variation

Realizing long-term monitoring of benthic habitat type variation has emerged as a crucial step to identify coral reef ecosystems at a risk of significant, detrimental transformation period, and to provide fundamental instruction for future ecosystem risk assessments and regulation (Murray et al. 2018). Spaceborne earth observation has become the most effective data source in observing the spatial distribution and temporal variation of reef ecosystems, particularly for those remotely and hard accessible located away from land (Foo and Asner 2019; Purkis et al. 2019). In this context, the ability of data processing framework delivered in this work to provide coincident geomorphic zonation map and benthic habitat map was of great advantage. As shown in Figure 8, according to Madin and Madin (2015), consolidated reefs are suitable basis for coral reef to brought up; thus, the predominant geomorphic zonation appropriate for coral growth would be reef slope, reef front, reef crest, and outer reef flat area, which is consistent with coral distribution results presented in this research. Among all benthic habitat categories, seagrass had the lowest classification accuracy, followed by microalgal mats. Compared to other benthic habitat type, the proportion of seagrass and microalgal mats was relatively small and the spectrum of them was easily confused with adjacent benthic habitat type conceiving the spatial resolution is 30 m of Landsat-8 OLI images. Moreover, for water depth deeper than reef front, variation of DII decreased as bottom reflectance attenuation increased with increasing water depth; hence, dark pixels around reef front were all supposed as consolidated material (coral or rock). Despite good agreement between these two data sets, it did not mean that Landsat-8 OLI images could be used to pinpoint coral reef, considering its relative coarse spatial resolution. Here, it was an indication that under different spatial scales, coral abundance was of significant relation in similar magnitude.

In particular, both geomorphic zonation and benthic habitat classification maps established in the anterior steps were the results in 2020. Considering the relative stability of geomorphology characteristics, geomorphic zonation map could be regarded as an overall state in a short time. Nevertheless, benthic habitat was susceptible to environmental changes; thus, the temporal variation of benthic types needs to be further explored. Initially, the proportion of each benthic type at each reef of the Xisha Islands in 2020 was depicted in Figure 11. Benthic distribution showed different features in different reefs. Limited by sensor’s spatial resolution, seagrass and microalgal mats were barely distinguished separately and generally mixed with other benthic types. Therefore, the classification of benthic habitat mainly focused on coral, sand, rock, and rubble that have distinct spectral or morphological differences. Since the benthic habitat classification results in 2020 were with acceptable accuracy, a natural question arose that whether the distribution of coral in different years (here we only took area of coral into consideration) retrieved based on the above procedure could accurately reflect the long-term status of coral and also hint the potential trend in the future? The clear Landsat-8 OLI images in different years at each reef location were considered to accomplish benthic classification step and eight reefs/islands were retained having continuous time series data. The area of coral was determined and compared with degree heat week (DHW) at reefs independently over a 9-year period (2013–2021) (Figure 12), which covered six El Niño and La Niña events in these regions.

Figure 11. Distribution map of benthic habitat components at each reef in the Xisha Islands in 2020. Percentage of land cover at each reef is presented with green bars.

Figure 12. Coral area (blue dots) at several reefs in the Xisha Islands and DHW (grey lines with black dots) from 2013 to 2021. El Niño are highlighted with color presented strength (weak=yellow; moderate=orange; strong=red), and La Niña are highlighted with color presented strength (weak=light blue; moderate=sky blue; strong=steel blue).

Maximum heat stress was detected in northwest part of the Xisha Island, where DHW reached around 13.5°C-week (Bei Reef, Xuande Reefs, and Jinyin Island), and intensity gradually decreased toward southeast direction, supporting previous findings by Lu et al. (2022). Common fluctuation in area of corals was found in most regions that area presented with steadily decreasing trend between 2013 and 2016, followed by a slow increase from 2016 to 2019 and a sharp decline after 2019, which overall was consistent with the evidence provided by Wu et al. (2022). From 2014 to 2015, El Niño events occurred and developed from weak to strong, and the El Niño-driven atmospheric and oceanic changes (e.g. cloud cover, humidity, monsoon wind) would regulate sea surface temperature (SST) of SCS (Wang et al. 2006). The SCS SST often increased under the influence of enhanced subsidence, reduced cloud cover, and increased solar radiation absorbed by the ocean during El Niño events, and the maximum SST anomaly generally appeared lagging the El Niño peak phase by 3–6 months (Klein, Soden, and Lau 1999). Therefore, two DHW peaks emerged in succession from 2014 to 2015, indicating two thermal stress incidents in the Xisha Islands, which may cause coral reef bleaching in this region. Another obvious DHW peak was discovered during 2020 accompanied by the La Niña incident. According to Chen et al. (2022), anomalous anticyclones over the north SCS associated with the Indian Ocean warming and La Niña aroused extreme warming of SCS and directly resulting in abnormal high DHW value. When DHW reached the maximum value, potential coral bleaching, even death, could be expected. Turning to the consequences of area variation during these two periods, coral area decreased slightly during 2014 and 2015, and several regions experienced a rapid rebound subsequently (e.g. Xuande Reefs, Jinyin Island, Yuzhuo Reef, and Panshiyu Island). However, almost all reef regions went through a drastical decline in 2020, and the subsequent restoration or more serious degradation was not clear yet limited by the study period. Comparing the performance of coral area changes in these two stages, it could be speculated that corals in the Xisha Islands have good resistance and tolerance to chronic heat stress in 2014–2015, when two successive relatively mild DHW peaks appeared, while its systemic resilience has limitation and may be abrupt in acute heat stress events (Cheung et al. 2021).

In this study, only sea surface temperature-based DHW was employed here as an environmental indicator, and area change of corals related with DHW implied that corals in the Xisha Islands may possess the ability to withstand chronic thermal stresses, yet still not capable of adapting to acute heat stress, which requires further evidence. Moreover, several other related environmental variables, such as pH, photosynthetically active radiance (PAR), turbidity, chlorophyll-a, fishing pressure, population density, crown of thorns starfish eruption, etc. (Hoegh-Guldberg 1999; Mason, Skirving, and Dove 2020) were not incorporated and discussed here. According to Galbraith and Convertino (2021), variation of coral area also could be seen as a signal emitted by the microbiome responding to biological organization and environmental pressure. Thus, the coral area variation found in this research to some extent may also be the result of multivariable compound effect. More attention should be paid in understanding subtle eco-environmental interaction patterns to offer more objective and reliable reference for coral reef monitoring and early warning. Under such consideration, eco-environmental factor predictive causality could be further explored based on spatial and biogeochemical network inference model, combined with sensitivity analysis methods (Wang, Galbraith, and Convertino 2023; Pianosi et al. 2016; Runge et al. 2019), which will not only help in pinpointing risk coral areas, but also determine environmental triggers and assist managers in targeted risk control and source prevention. Overall, coral area fluctuation derived from benthic habitat classification results provided in this research reflected that coral reef in the Xisha Island may own the ability to resist chronic heat stress, while extra attention and timely protective measures are required when acute heat stress incidents occur.

4.4 Limitations and future work

The proposed framework in this study poses a few restrictions. First of all, the accuracy of PI results was confined with data acquisition frequency, visible light penetration depth, and high instability geomorphology of some sands. Second, we applied 30 m resolution satellite images to determine geomorphic zonation and benthic habitat distribution, while remote sensing images of high spatial resolutions would be more appropriate in comprehensive understanding of coral reefs at scales of a few meters or even less. Third, the results were compared and validated with historical record data and limited measured data. Last but not least, the driving mechanism of coral area variation in the Xisha Islands is complex; therefore, the inter actions between the coral areas and biogeochemical factors (ecosystem connectome) should not be underestimated. In future work, the combination of more in-situ data, Landsat series data, and other remote sensing images with higher spatial and temporal resolutions (e.g. Sentinel-2) hopefully can provide improved local to global mapping and monitoring of coral reef ecosystems at finer scale. Moreover, it is also essential to consider multiple stressors to coral reef health and interaction patterns between coral reefs and its surrounding environment; thus, dynamic models or ecological assessment model would be considered as a useful tool in future research to assess the synergistic effects of local and global stressors on coral reef ecosystem functions.

5. Conclusion

Coral reefs are increasingly endangered under multiple disturbances from both changing environment and anthropogenic interference. In order to maintain the coral reef ecosystem and strengthen their resilience to external changes, a detailed exploration of reef state, like reef extent, geomorphic zonation, and benthic habitat mapping of coral reef is fundamental and essential. This research established an agile but effective processing framework accomplishing rapid geomorphologic and benthic habitat types classification of the Xisha Islands with unsupervised classification algorithms based on a combination of time-series Landsat-8 satellite data without in-situ measurements. A concrete and original inventory including reef spatial extent, PI map, potential geomorphic zonation map, and benthic habitat map of the Xisha Islands (time averaged results) has been set up. In contrast to open-access datasets (UNEP-WCMC/RESDC), reef extent extracted in this research presented with better authenticity and boundary alignment when compared with RGB images from Landsat-8 in different years. Additionally, validation of geomorphic zonation/benthic habitat classification results (accuracy reached 72% and 86%, respectively) suggested that the outputs of this research are appropriate translation of remote sensing data from producers to end-users (reef managers, policy makers, conservation practitioners), which could be served as a reference in coral reef protection, regulation, management, and restoration work. Interestingly, the ability to resist chronic heat stress was found in the coral reef at Xisha Islands as the coral area presented a growth trend after two successive heat stress events in 2014–2015, while the potential recovery situation after acute heat stress incident occurred in 2020 still needs further follow-up studies. This study presents an attempt to rapidly map the current comprehensive situation of coral reef at Xisha Islands based only on publicly available remote sensing data with minimal processing limitations, and accuracy of results was acceptable under the premise of lack of measured data. We anticipate that the maps provided with this study could support specific localized uses, and integrated with multiple resources to enable wall-to-wall coverage for current local to global-scale sustainable management and preservation frameworks.

Acknowledgement

The authors would like to express sincere gratitude to Hainan Provincial Observatory of Ecological Environment and Fishery Resource in Yazhou Bay for providing field data, and the United States Geological Survey (USGS) for Landsat-8 data contributing. The authors are also grateful to the Allen Coral Atlas for coral reef products sharing, which help us lot in results comparison in this work.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/15481603.2023.2261213

Data availability statement

The data that support the findings of this study are available from the corresponding author, S.H., upon reasonable request.

Additional information

Funding

The research was supported by the Hainan Provincial Joint Project of Sanya Yazhou Bay Science and Technology City (Grant No. 2021JJLH0053), the Hainan Provincial Joint Project of Sanya Yazhou Bay Science and Technology City (Grant No. 2021CXLH0020), the Research Startup Funding from Hainan Institute of Zhejiang University (Grant No. 0208-6602-A12202).