Code and data: Spatial variation in upper limits of coral cover on the Great Barrier Reef

Identifying the maximum coral cover that a coral community can sustain (i.e., its ‘upper limit’) is important for predicting community dynamics and improving management strategies. Here, we quantify the relationship between estimated upper limits and key environmental factors on coral reefs: hard substrate availability, temperature, and water clarity. We used 32 years (1990-2022) of data on coral cover around reef perimeters in the Great Barrier Reef, Australia. Each reef was divided into four wave exposure habitats depending on prevailing wind conditions. For each site, we determined if hard coral cover had reached a plateau or upper limit. Next, we extracted existing estimates of hard substrate availability, modelled water temperature, and Secchi depth. Then, we quantified the relationship between these environmental variables and the upper limits. We found varying upper limits across the GBR, with a median of 33% coral cover and only 17% of the estimated upper limits exceeded 50% coral cover. Upper limits increased towards the southern reefs. Our results show that upper limits increased with increasing hard substrate availability and decreased with temperature and, to a lesser extent, with water clarity. The upper limits estimated in this study are much lower than what is commonly assumed when modelling ecological dynamics, most likely resulting in predicted recovery rates being inappropriately high. Although hard substrate ultimately restricted upper limits, there are mechanisms constraining the proportion of hard substrate that is covered by hard corals. The negative relationship between temperature and upper limits cannot be explained by changes in macroalgal abundance but may be related to changes in species composition. The quantitative relationships between the upper limits of coral cover and environmental variables will provide critical information to prioritise sites for management interventions.

Methods
Coral cover data

We used hard coral cover data from manta tow surveys conducted by the Australian Institute of Marine Science’s (AIMS) Long-term Monitoring Program (LTMP) from 1990 to 2022 (Emslie et al., 2020). Manta tow surveys quantified broadscale patterns in percent coral cover on a ten-meter swath of the reef slope around the entire perimeters of 173 reefs on the Great Barrier Reef, Australia (from 12°S to 24°S and 143°E to 152°E). Reefs were visited a maximum of once a year.

Each reef was divided into four reef zones or sites based on prevailing wind exposure: one back (leeward) site that was the most protected, one front (windward) site being the most exposed, and two flank sites with intermediate exposure. Sites were classified based on reef geomorphology and local knowledge of prevailing wind conditions. Each tow was assigned a reef zone, and from 2014 onwards, the tow paths were georeferenced. Depending on weather conditions and variable reef structure, following the exact path was not always possible. Since community composition changes with wave exposure, we looked for upper limits in each of the four sites within reef separately. 

Observers were towed by a boat at constant speed, and for each two-minute tow, visually estimated the proportion of the reef slope covered by live hard corals using a categorical coral cover score (0, >0-5%, >5-10%, >10-20%, >20-30%, >30-40%, >40-50%, >50-62.5%, >62.5-75%, >75-87.5%, >87.5-100%) (Miller et al., 2018). During the two-minute tow, the observer is towed parallel to the reef crest, aiming to visually scan a ten-meter band of the reef slope, just below the reef crest. The length of the transect varied depending on the distance covered during the two-minute duration of the tow, with a mean tow length of 200 m (with 95% of the tows between 134 – 269 m). Only one categorical hard coral cover score was recorded for each two-minute tow.

For analyses, coral cover for each two-minute tow was converted to the midpoint of the cover category recorded. For each year surveyed, all tows from each site within reef were summarised by calculating the median to minimise the effect of outliers. A small subset of reefs had inconsistent wave exposure zone assignations that could not be corrected using the georeferenced tows and therefore we excluded sites within reefs that were problematic (three entire reefs were excluded, and 13 additional reefs had at least one site excluded).

To estimate the upper limit of coral cover for each site within each reef, we first calculated the five-year running mean of the median coral cover. From the five-year running mean, we identified instances when coral cover reached a plateau. We called coral cover at the plateau the ‘upper limit’, and we defined it as the maximum running mean cover over a period of five undisturbed years. The upper limit had to: 1) be within 80% of the maximum running mean for that site across all surveyed years, 2) the two previous and the two following years also had a running mean within 80% of the maximum value, and 3) the running mean estimate was calculated with at least three data points (i.e., only two years out the five could have an interpolated value). We chose these criteria to ensure that we only captured instances of stable high coral cover within each site, but with some flexibility (e.g., in terms of interpolation) so that we had enough upper limit estimates to capture their variation across space. Since the specific criteria are somewhat arbitrary, we repeated the process of estimating upper limits with different running mean windows (5, 7, 9, and 11 years) and different thresholds for the difference between the running mean and the maximum coral cover (80, 85, and 90%). In total, we had 12 sets of upper limit estimates, each estimated using different criteria. For any definition of upper limits, each reef could have a maximum of four upper limit values, one for each site (one exposed, one sheltered, and two intermediate reef sites). 

Environmental variables

We extracted long-term medians of the environmental values for hard substrate availability, temperature, and water clarity, assuming it to be representative of the general environmental conditions across sites. 

We extracted estimates of the proportion of reef surfaces constituting hard substrate, derived from freely available benthic maps known as “Great Barrier Reef 10m Grid GBRMP Benthic” (Lyons et al., 2020; Great Barrier Reef Marine Park Authority, 2021). The maps are generated using high spatial resolution satellite imagery, bathymetry, wave height estimates, and field data for training and validation (Roelfsema et al., 2021). In the maps, each 10m-by-10m of planar area pixel down to 10m mean sea level depth is assigned one of four types of substrate: coral or macroalgae, rock, rubble, or sand. Coral cannot be distinguished from macroalgae in the satellite imagery, and any pixels labelled as coral/macroalgae or rock (i.e., hard substrate) were considered suitable habitat for corals to establish and grow, unlike rubble and sand. 

Each reef site had multiple pixels with a benthic composition label. To ensure we captured the benthic composition of the same locations that were surveyed with the manta tows, we divided each reef into 100m x 100m squared grid cells, and we assigned each cell with a georeferenced tow to a reef site. For each reef site, we extracted benthic composition estimates from all assigned cells and determined the proportion of 10m x 10m pixels labelled as hard substrate divided by the total number of pixels.

To investigate the relationship between upper limits and key environmental variables, we summarised long-term medians in temperature and Secchi depth (a measure of water clarity, with values above 10 - 15 m representing high water clarity) from hydrodynamic and biogeochemical models. Daily sea water temperatures at a depth of 7.5m below mean sea level and Secchi depth were extracted from the eReefs marine models (GBR1 v2.0 and GBR4, respectively) (Steven et al., 2019) using the R-package “ereefs” (Robson, 2023). eReefs is an information system that predicts biogeochemical and hydrodynamic variables across the Great Barrier Reef. Daily estimates were available at a spatial resolution coarser than the reef site (~1 x 1 km and ~4 x 4 km) during a period of six and eight years for temperature and Secchi depth, respectively (temperature: 1 January 2015 to 1 January 2021, Secchi depth: 1 December 2010 to 1 December 2018). For these variables, we used the georeferenced tows to find the latitude and longitude at the centre of the reef site and computed the median of the daily estimates for the grid cell upon which the centre of each reef site fell. The eReefs hydrodynamic models consist of a nominal 1 km resolution hydrodynamic model covering the Great Barrier Reef World Heritage Area nested within a nominal 4 km resolution regional model that extends across the Coral Sea, which is in turn nested within a global circulation model. Among its outputs, the hydrodynamic model generates predictions for water temperature using the three-dimensional, baroclinic, finite difference hydrodynamic model Sparse Hydrodynamic Ocean Code (SHOC; Herzfeld, 2006).  Although reefs differed in their depth distribution and geomorphology, we expect a depth of 7.5m to overlap with all manta tows swaths. Estimates of Secchi depth were obtained from the eReefs biogeochemical model that has a 4 km resolution. The optical model used to generate Secchi depth predictions can reproduce remote-sensing reflectance patterns that arise from specific biogeochemical states (Baird et al., 2016). Secchi depth was included as a proxy for water clarity and net productivity.

Statistical analyses

All analyses were performed in R version 4.3.1 (R Core Team, 2023). We fitted a series of linear mixed effects models in a Bayesian framework using the R package “brms” version 2.20.4 (Bürkner, 2017) to quantify how sites’ upper limit estimates related to environmental conditions. Because reefs could have up to four data points each (in different wave exposure sites), we included a random effect of reef in the models. Because our response variable - upper limit – was a proportion (converted to percentage in the figures), we fitted the models assuming a beta distribution.

For each of the 12 sets of upper limit estimates (i.e., one for each running mean window and percentage threshold combination), we fitted seven models to test whether the combination of fixed effects affected the magnitude and direction of the environmental variable coefficients. There was one model that included all three environmental variables (the proportion of hard substrate, temperature, and Secchi depth) as fixed effects, three models that included two out of the three environmental variables only, and three models that included only one of the environmental variables. This gave a total of 84 fitted models.

In each model, the environmental variables were standardised as Z-scores ( , where  is the value of the environmental variable for that site, and  and  are the mean and standard deviation across all observations, respectively). For simplicity and since one third of the sets of upper limit estimates had less than 50 estimates, we avoided fitting interactions between environmental variables.

Additionally, to test whether the relationship between coral cover and the three environmental variables differed among high, medium, and low coral cover at a site, we fitted three quantile regressions. Each quantile regression represented high (95th quantile), medium (50th quantile), and low (10th quantile) coral cover. The response variable was the yearly coral cover at each site (n = 7787). Temperature, hard substrate availability, and Secchi depth (as Z-scores) were included as the fixed explanatory variables and, since there were multiple data points per reef site, site nested within reef was included as a random effect.

For all models, the priors were specified by using the function brms ‘get_prior’, which gave a flat prior for the fixed effects and intercept and weakly informative priors for the hierarchical parameters (a Student’s t-distribution with a normality parameter of 3, with a mean of 0, and a standard deviation of 2.5 for the models with a beta distribution and a standard deviation of 14.8 for the quantile regressions). Each model had three chains of 20,000 iterations (half of which were discarded as warm up) and a thinning of 5, resulting in 6,000 posterior distribution draws. To investigate how consistent was the relationship between each environmental variable and coral cover upper limits, we recorded the coefficient estimate for each variable included in each model and whether its 95% credible interval overlapped with zero.

Please refer to the published manuscript for more detailed information on data and processing.

References

Bürkner, P.-C. (2017) brms: An R Package for Bayesian Multilevel Models Using Stan . Journal of Statistical Software , 80, 1–28.

Emslie, M.J., Bray, P., Cheal, A.J., Johns, K.A., Osborne, K., Sinclair-Taylor, T. & Thompson, C.A. (2020) Decades of monitoring have informed the stewardship and ecological understanding of Australia’s Great Barrier Reef. Biological conservation, 252, 108854.

Great Barrier Reef Marine Park Authority (2021) GBR10 GBRMP Benthic.

Herzfeld, M. (2006) An alternative coordinate system for solving finite difference ocean models. Ocean Modelling, 14, 174–196.

Lyons, M.B., Roelfsema, C.M., Kennedy, E. V, Kovacs, E.M., Borrego‐Acevedo, R., Markey, K., Roe, M., Yuwono, D.M., Harris, D.L. & Phinn, S.R. (2020) Mapping the world’s coral reefs using a global multiscale earth observation framework. Remote Sensing in Ecology and Conservation, 6, 557–568.

Miller, I.R., Jonker, M.J. & Coleman, G. (2018) Crown-of-thorns starfish and coral surveys using the manta tow technique, Townsville, Australia.

R Core Team (2023) R: A language and environment for statistical computing.

Robson, B. (2023) ereefs: Useful Functions to Handle eReefs and EMS model Output.

Roelfsema, C.M., Lyons, M.B., Castro-Sanguino, C., Kovacs, E.M., Callaghan, D., Wettle, M., Markey, K., Borrego-Acevedo, R., Tudman, P. & Roe, M. (2021) How Much Shallow Coral Habitat Is There on the Great Barrier Reef? Remote Sensing, 13, 4343.

Steven, A.D.L., Baird, M.E., Brinkman, R., Car, N.J., Cox, S.J., Herzfeld, M., Hodge, J., Jones, E., King, E. & Margvelashvili, N. (2019) eReefs: An operational information system for managing the Great Barrier Reef. Journal of Operational Oceanography, 12, S12–S28.

研究背景
明确珊瑚群落所能维持的最大珊瑚覆盖度（即其‘上限’），对于预测群落动态、优化管理策略具有重要意义。本研究量化了珊瑚礁上估算的覆盖度上限与关键环境因子之间的关系：硬质基底可用性、海水温度与水体透明度。我们使用了澳大利亚大堡礁（Great Barrier Reef, GBR）1990-2022年共32年的礁周珊瑚覆盖度数据。依据盛行风条件，每个礁体被划分为4个波浪暴露生境。针对每个样点，我们判断其硬质珊瑚覆盖度是否已达到平台期或上限。随后，我们提取了现有硬质基底可用性、模拟海水温度及赛氏盘深度（Secchi depth）的估算值，并量化了这些环境变量与覆盖度上限之间的关系。研究发现，大堡礁范围内的覆盖度上限存在显著空间异质性，中位值为33%的珊瑚覆盖度，仅17%的估算上限超过50%珊瑚覆盖度，且上限值随礁体向南偏移而升高。结果表明，覆盖度上限随硬质基底可用性提升而升高，随海水温度升高、以及一定程度上随水体透明度下降而降低。本研究估算的覆盖度上限远低于生态动力学建模中的常规假设，这大概率会导致预测的珊瑚恢复速率过高。尽管硬质基底最终限制了覆盖度上限，但仍存在多种机制约束着硬质基底被硬质珊瑚占据的比例。海水温度与覆盖度上限之间的负相关关系无法通过大型藻类丰度的变化解释，但可能与珊瑚物种组成的改变有关。珊瑚覆盖度上限与环境变量之间的定量关系，将为优先确定管理干预样点提供关键依据。

方法
1. 珊瑚覆盖数据
我们使用了澳大利亚海洋科学研究所（Australian Institute of Marine Science, AIMS）长期监测计划（Long-term Monitoring Program, LTMP）于1990-2022年开展的蝠鲼拖曳调查（manta tow surveys）所获取的硬质珊瑚覆盖度数据（Emslie等，2020）。该调查针对澳大利亚大堡礁173个礁体（纬度范围12°S至24°S，经度范围143°E至152°E）的整个礁周斜坡的10米宽条带，量化了珊瑚覆盖度的大尺度格局。每个礁体每年最多被调查一次。
依据盛行风暴露程度，每个礁体被划分为4个礁区或样点：1个背风（庇护）样点为受保护程度最高的区域，1个迎风（暴露）样点为受波浪冲击最强的区域，以及2个具有中等暴露程度的侧翼样点。样点划分基于礁体地貌与盛行风的本地认知。每个拖曳采样点被分配至对应的礁区，2014年起，拖曳路径被进行地理配准。受天气条件与复杂礁体结构的影响，严格遵循预设路径并不总能实现。由于群落组成随波浪暴露程度变化，我们分别在礁体内的4个样点中分别寻找覆盖度上限。
调查人员以恒定速度乘船拖曳采样，每次2分钟的拖曳过程中，通过分类珊瑚覆盖度评分（0、>0-5%、>5-10%、>10-20%、>20-30%、>30-40%、>40-50%、>50-62.5%、>62.5-75%、>75-87.5%、>87.5-100%）目视估算礁坡被活硬质珊瑚覆盖的比例（Miller等，2018）。2分钟拖曳过程中，调查人员沿礁脊平行航行，旨在目视扫描礁脊正下方10米宽的礁坡条带。拖曳的长度随2分钟航行过程中覆盖的距离变化，平均拖曳长度为200米（95%的拖曳长度介于134-269米之间）。每次2分钟的拖曳仅记录一个分类硬质珊瑚覆盖度评分。
为便于分析，我们将每次2分钟拖曳的珊瑚覆盖度转换为对应覆盖度分类的中点值。对于每一个调查年份，我们通过计算中位数来汇总礁体内每个样点的所有拖曳数据，以降低异常值的影响。一小部分礁体的波浪暴露区划分存在不一致，且无法通过地理配准的拖曳数据修正，因此我们排除了存在问题的礁体内的样点：共排除3个完整礁体，另有13个礁体各排除至少1个样点。
为估算每个礁体内每个样点的珊瑚覆盖度上限，我们首先计算了中位数珊瑚覆盖度的5年滑动平均值。从5年滑动平均序列中，我们识别出珊瑚覆盖度达到平台期的时段，将该平台期的珊瑚覆盖度定义为‘上限’，即5个未受扰动年份内的最大滑动平均覆盖度。该上限需满足以下3个条件：1）在该样点所有调查年份的最大滑动平均覆盖度的80%范围内；2）前后各2年的滑动平均覆盖度也处于最大值的80%范围内；3）滑动平均估算值至少基于3个数据点（即5年中仅可有2年为插值得到的值）。我们设置这些标准以确保仅捕获每个样点内稳定的高珊瑚覆盖度时段，同时保留一定灵活性（例如插值规则），从而获得足够的上限估算值以反映空间变异。由于具体标准存在一定主观性，我们采用不同的滑动平均窗口（5、7、9、11年）以及滑动平均与最大珊瑚覆盖度的差异阈值（80%、85%、90%）重复了上限估算流程，最终得到12组上限估算结果，每组均基于不同的标准。对于任意一种上限定义方式，每个礁体最多可拥有4个上限值，对应4个样点（1个暴露样点、1个庇护样点与2个中等暴露样点）。

2. 环境变量
我们提取了硬质基底可用性、温度与水体透明度的长期中位数，以此代表样点的一般环境条件。
我们提取了硬质基底占比的估算值，该数据来自免费公开的底栖地图——‘大堡礁10米网格GBRMP底栖地图’（Great Barrier Reef 10m Grid GBRMP Benthic）（Lyons等，2020；大堡礁海洋公园管理局，2021）。该地图通过高空间分辨率卫星影像、水深数据、波浪高度估算值与野外调查数据进行训练与验证生成（Roelfsema等，2021）。在地图中，平均海平面以下10米深度范围内的每一个10m×10m平面像素被划分为4类基底：珊瑚或大型藻类、岩石、碎石或沙质。由于卫星影像无法区分珊瑚与大型藻类，所有被标记为珊瑚/大型藻类或岩石的像素（即硬质基底）均被视为珊瑚定植与生长的适宜生境，而碎石与沙质则不属于此类。
每个礁体样点拥有多个带有底栖组成标签的像素。为确保我们获取的底栖组成与蝠鲼拖曳调查的采样位置一致，我们将每个礁体划分为100m×100m的网格单元，并将带有地理配准的拖曳采样点分配至对应的礁区网格单元。针对每个礁体样点，我们从所有分配的网格单元中提取底栖组成估算值，并计算标记为硬质基底的10m×10m像素占总像素数的比例。
为探究覆盖度上限与关键环境变量之间的关系，我们从水动力与生物地球化学模型中提取了温度与赛氏盘深度（水体透明度的衡量指标，10-15米以上代表高水体透明度）的长期中位数。我们通过R包‘ereefs’（Robson，2023）从eReefs海洋模型（分别为GBR1 v2.0与GBR4）中提取了平均海平面以下7.5m深度处的每日海水温度与赛氏盘深度数据（Steven等，2019）。eReefs是一套可预测大堡礁范围内生物地球化学与水动力变量的信息系统。每日估算数据的空间分辨率略低于礁体样点（温度数据约为1×1km，赛氏盘深度数据约为4×4km），数据覆盖时段分别为：温度：2015年1月1日至2021年1月1日；赛氏盘深度：2010年12月1日至2018年12月1日。针对这些变量，我们通过地理配准的拖曳采样点获取礁体样点中心的经纬度，并计算该中心所在网格单元内每日估算值的中位数。eReefs水动力模型采用标称1km分辨率的水动力模型嵌套于覆盖珊瑚海的标称4km分辨率区域模型中，而该区域模型又嵌套于全球环流模型中，该水动力模型覆盖大堡礁世界遗产区域。该模型的输出之一为通过三维斜压有限差分水动力模型Sparse Hydrodynamic Ocean Code（SHOC; Herzfeld, 2006）生成的海水温度预测值。尽管不同礁体的深度分布与地貌存在差异，但我们预期7.5m深度可覆盖所有蝠鲼拖曳采样的礁坡条带。赛氏盘深度估算值来自分辨率为4km的eReefs生物地球化学模型，用于生成赛氏盘深度预测值的光学模型可复现由特定生物地球化学状态引起的遥感反射率模式（Baird等，2016）。赛氏盘深度被用作水体透明度与净生产力的替代指标。

3. 统计分析
所有分析均在R 4.3.1版本中完成（R核心团队，2023）。我们使用R包‘brms’2.20.4版本（Bürkner, 2017）在贝叶斯框架（Bayesian framework）下拟合了一系列线性混合效应模型（linear mixed effects models），以量化样点的覆盖度上限估算值与环境条件之间的关系。由于每个礁体最多可拥有4个数据点（对应不同波浪暴露样点），我们在模型中加入了礁体的随机效应。由于响应变量——覆盖度上限——为比例值（在图中转换为百分比），我们采用Beta分布（Beta distribution）进行模型拟合。
针对12组上限估算结果（即每组对应不同的滑动平均窗口与百分比阈值组合），我们分别拟合了7个模型，以检验固定效应组合如何影响环境变量系数的大小与方向。其中1个模型包含全部3个环境变量（硬质基底占比、温度与赛氏盘深度）作为固定效应，3个模型仅包含3个变量中的2个，另有3个模型仅包含1个环境变量，总计拟合84个模型。
在每个模型中，环境变量均被标准化为Z分数（Z-scores，公式为$Z = frac{x-mu}{sigma}$，其中$x$为该样点的环境变量值，$mu$与$sigma$分别为所有观测值的均值与标准差）。为简化分析，且由于三分之一的上限估算组的有效估算值少于50个，我们未加入环境变量间的交互项。
此外，为检验珊瑚覆盖度与3个环境变量之间的关系是否在样点的高、中、低覆盖度水平上存在差异，我们拟合了3个分位数回归（quantile regression），分别代表高覆盖度（95%分位数）、中等覆盖度（50%分位数）与低覆盖度（10%分位数）。响应变量为每个样点的年度珊瑚覆盖度（样本量n=7787）。固定解释变量为标准化为Z分数的温度、硬质基底可用性与赛氏盘深度，由于每个礁体样点拥有多个数据点，我们加入了‘样点嵌套于礁体’的随机效应。
对于所有模型，我们通过brms的‘get_prior’函数设置先验分布：固定效应与截距项采用平坦先验，层级参数采用弱信息先验（对于Beta分布模型，采用自由度为3的学生t分布，均值为0，标准差为2.5；对于分位数回归模型，标准差为14.8）。每个模型运行3个马尔可夫链，每条链包含20000次迭代（其中一半作为预热期被舍弃），并以5为间隔进行thinning，最终得到6000个后验分布抽样结果。为探究每个环境变量与珊瑚覆盖度上限之间的关系是否具有一致性，我们记录了每个模型中所包含变量的系数估算值，以及其95%可信区间是否与0存在重叠。
如需获取数据与处理流程的更多详细信息，请参阅已发表的论文原稿。

参考文献
Bürkner, P.-C. (2017) brms：一款基于Stan的贝叶斯多层模型R包. 《统计软件杂志》, 80, 1–28.
Emslie, M.J., Bray, P., Cheal, A.J., Johns, K.A., Osborne, K., Sinclair-Taylor, T. & Thompson, C.A. (2020) 数十年监测助力澳大利亚大堡礁的管理与生态认知. 《生物保护》, 252, 108854.
大堡礁海洋公园管理局 (2021) GBR10 GBRMP底栖地图.
Herzfeld, M. (2006) 有限差分海洋模型求解的替代坐标系. 《海洋建模》, 14, 174–196.
Lyons, M.B., Roelfsema, C.M., Kennedy, E. V, Kovacs, E.M., Borrego‐Acevedo, R., Markey, K., Roe, M., Yuwono, D.M., Harris, D.L. & Phinn, S.R. (2020) 基于全球多尺度地球观测框架绘制全球珊瑚礁地图. 《遥感在生态学与保护中的应用》, 6, 557–568.
Miller, I.R., Jonker, M.J. & Coleman, G. (2018) 采用蝠鲼拖曳技术开展的长棘海星与珊瑚调查, 澳大利亚汤斯维尔.
R核心团队 (2023) R：统计计算的语言与环境.
Robson, B. (2023) ereefs：用于处理eReefs与EMS模型输出的实用函数.
Roelfsema, C.M., Lyons, M.B., Castro-Sanguino, C., Kovacs, E.M., Callaghan, D., Wettle, M., Markey, K., Borrego-Acevedo, R., Tudman, P. & Roe, M. (2021) 大堡礁的浅海珊瑚栖息地面积有多大？《遥感》, 13, 4343.
Steven, A.D.L., Baird, M.E., Brinkman, R., Car, N.J., Cox, S.J., Herzfeld, M., Hodge, J., Jones, E., King, E. & Margvelashvili, N. (2019) eReefs：用于大堡礁管理的业务化信息系统. 《业务海洋学杂志》, 12, S12–S28.