Mathematical modelling and uncertainty quantification for analysis of biphasic coral reef recovery patterns

Abstract

Coral reefs are increasingly subjected to major disturbances threatening the health of marine ecosystems. Substantial research is underway to develop intervention strategies that assist reefs in recovery from, and resistance to, inevitable future climate and weather extremes. To assess potential benefits of interventions, mechanistic understanding of coral reef recovery and resistance patterns is essential. Recent evidence suggests that more than half of the reefs surveyed across the Great Barrier Reef (GBR) exhibit deviations from standard recovery modelling assumptions when the initial coral cover is low ($\leq 10$%). New modelling is necessary to account for these observed patterns to better inform management strategies. We consider a new model for reef recovery at the coral cover scale that accounts for biphasic recovery patterns. The model is based on a multispecies Richards’ growth model that includes a change point in the recovery patterns. Bayesian inference is applied for uncertainty quantification of key parameters for assessing reef health and recovery patterns. This analysis is applied to benthic survey data from the Australian Institute of Marine Science (AIMS). We demonstrate agreement between model predictions and data across every recorded recovery trajectory with at least two years of observations following disturbance events occurring between 1992–2020. This new approach will enable new insights into the biological, ecological and environmental factors that contribute to the duration and severity of biphasic coral recovery patterns across the GBR. These new insights will help to inform managements and monitoring practice to mitigate the impacts of climate change on coral reefs.

Keywords: coral reef recovery; biphasic recovery patterns; reef recovery modelling; climate change; marine monitoring.

1 Introduction

Environmental stressors, including climate change, continue to drive increased frequency and severity of disturbance events affecting coral reefs worldwide. This is leading to long-term declines in coral cover and changes in community composition (De’ath et al.,, 2009, 2012; Duran et al.,, 2017; Johns et al.,, 2014; Lapointe et al.,, 2019; Tebbett et al.,, 2023). To ensure the future health of coral reefs and connected marine ecosystems we need to continue developing capabilities to effectively prioritise targeted intervention and restoration efforts (McLeod et al.,, 2022). Forecasts of coral reef health under different potential environmental scenarios are an essential part of this decision-making process to ensure limited resources are efficiently allocated (Thompson and Dolman,, 2010; Vercelloni et al.,, 2020; van Woesik et al.,, 2018). In addition, predictions of expected versus observed coral cover is an important health indicator of recovery  (Castro-Sanguino et al.,, 2021; Logan et al.,, 2020; Thompson et al.,, 2020). To produce reliable forecasts and predictions for decision-makers, we need models that are informed by our understanding of coral recovery patterns of reefs following perturbations by disturbance events, such as storm swells (Cheal et al.,, 2017; De’ath et al.,, 2012; Fabricius et al.,, 2008), predation by crown-of-thorns starfish, Acanthaster spp. (COTS) (De’ath et al.,, 2012; Pratchett et al.,, 2021; Seymour and Bradbury,, 1999; Vercelloni et al.,, 2017), and acute periods of thermal stress resulting in coral bleaching (Hughes et al.,, 2017; Mumby et al.,, 2021).

Recently, Warne et al., (2022) demonstrated evidence for biphasic recovery patterns amongst almost $50\%$ of reefs recovering from very low cover ($\leq 10\%$), across Australia’s Great Barrier Reef (GBR). Such recovery patterns were characterised by a transient period of reduced recovery that persists for approximately four years on average (Warne et al.,, 2022). Such recovery patterns represent deviations from standard recovery modelling approaches that are based on variations of classical population growth models, such as exponential, Gompertz, or Logistic growth models (MacNeil et al.,, 2019; Ortiz et al.,, 2018; Simpson et al.,, 2022; Thompson and Dolman,, 2010; van Woesik et al.,, 2018). Warne et al., (2022) also demonstrate coral recovery forecasts that ignore biphasic patterns could be overestimating future median cover estimates by up to 22% in absolute coral cover. Since the statistical methodology applied by Warne et al., (2022) is based upon multi-segment regression analysis (Aminikhanghahi and Cook,, 2016; Hinkley,, 1970; Jin et al.,, 2017), the study was restricted to a small subset of severely disturbed reefs with at least 5 observations of recovery. This motivates the development of a recovery modelling framework that can account for the effects of biphasic recovery in a broader range of scenarios.

In the biological and ecological modelling literature, biphasic models have been widely studied and applied in a variety of settings (Bodnar et al.,, 2013; Brouwer et al.,, 2017; Honsey et al.,, 2016; de A. Møller et al.,, 2021; Phaiboun et al.,, 2015), the most general of which has been demonstrated recently by Murphy et al., (2022). Here, we present a multi-species biphasic coral recovery model building on the modelling frameworks of Murphy et al., (2022) and Simpson et al., (2023). Using Bayesian analysis for model calibration, prediction, and uncertainty quantification, we demonstrate the approach is directly applicable to inform health assessment of the reefs within the GBR based upon coral cover monitoring data (Jonker et al.,, 2008). Through embedding our recovery model within a Bayesian analysis framework, we provide a powerful tool for uncertainty quantification of key recovery parameters and coral cover predictions that account for this parameter uncertainty.

2 Materials and methods

In this section we briefly describe the benthic survey data used in this work. A description of our mathematical model that captures biphasic recovery patterns is then presented along with the Bayesian methods for model calibration and assessment of model fitness.

2.1 Benthic survey data

Coral cover time-series data are obtained from benthic surveys undertaken by the Australian Institute of Marine Sciences (AIMS) Long Term Monitoring Program (LTMP) and Marine Monitoring Program (MMP) between 1992 to 2020. Combined, these collection programs include taxonomically classified cover data obtained from over 373 sites across 135 reefs. These monitored reefs are characteristic of the inshore, middle and outer shelf areas of the GBR (Figure 1(a)). The dataset includes a mixture of annual and biennial surveys at depths between 2m-9m. At each site, five permanent transects (dimension 5m$\times$50m for LTMP and 5m$\times$20m for MMP) are photographed at regular intervals (1m intervals for LTMP and 0.5m intervals for MMP). Coral cover is estimated and classified using a 5-point stencil overlaid in each image (See Jonker et al., (2008) for details on the data collection protocol). Example coral cover time series are shown for Agincourt 1 (Figure 1(b)), Thetford (Figure 1(c)) and Lady Musgrave reefs (Figure 1(d)).

Figure 1: (a) Location of reef sites that are monitored as part of the Long-Term Monitoring Program (LTMP) and the Marine Monitoring Program (MMP). The MMP monitors the inshore reefs (red), and the LTMP monitors the Middle (green) and Outer (blue) reef shelves. Example coral cover time-series are shown for (b) Agincourt 1, (c) Thetford, and (d) Lady Musgrave reefs with their locations indicated with arrows (black) in (A). Here, total hard coral cover (black solid), the fast growing Acropridae family (black dotted), and other hard corals (black dashed) are shown over time along with official disturbance event records (red dashed). Error bars indicate the standard error of the coral cover record based on the variation at the transect level.

Following Ortiz et al., (2018), and Warne et al., (2022), disturbances were identified by statistically significant decline ($p\leq 0.05$) in cover according to a one-sided paired t-test at the transect level ($n=5$) for consecutive visits (Figure 1B–D). A recovery trajectory is defined as the sequence of monitoring surveys between consecutive disturbances inclusive of the disturbance event defining the start of the recovery. Across the GBR, the patterns of coral recovery between disturbance events varies widely, as shown in the examples in Figure 1. This motivates a modelling approach that can account for deviations from standard growth patterns typically defined using a Gompertz or logistic model (Thompson et al.,, 2016, 2020; Simpson et al.,, 2023). One limitation of the methodology developed in Warne et al., (2022) was the restriction to recovery trajectories with low initial cover ($\leq 10\%$) and have $n\geq 5$ benthic survey visits before the next disturbance. Our new modelling framework reduces these restrictions, requiring only that the number of visits between subsequent disturbances is $n\geq 3$ with no restrictions on the initial coral cover are necessary. From the LTMP and MMP benthic surveys, we obtain $N=737$ recovery trajectories among $336$ reef sites and $124$ reefs across the GBR.

Throughout this manuscript, we denote $\mathbf{X}^{(i)}(t_{j})$ as the vector of taxonomically classified coral cover data at the reef site level for the $i$-th recovery trajectory, $i=1,2,\ldots,N$, as surveyed at time $t_{j}$, and $j=0,1,\ldots,n^{(i)}$. Here, $j=0$ relates to the initial visit in which the reef site was perturbed by a disturbance event and $n^{(i)}$ is the number of subsequent visits in the $i$-th recovery trajectory. Thus the $i$-th observed recovery trajectory is $\mathcal{D}^{(i)}=[\mathbf{X}^{(i)}(t_{0}),\mathbf{X}^{(i)}(t_{1}),\ldots,\mathbf{X}^{(i)}(t_{n^{(i)}})]$ containing the time-series of benthic survey data for that recovery trajectory.

2.2 Model development

Our mathematical model of coral reef recovery accounts for biphasic recovery patterns as identified in previous work (Warne et al.,, 2022), and allows for more general single-phase recovery patterns beyond standard approaches such as logistic or Gompertz growth models (Simpson et al.,, 2022; Thompson and Dolman,, 2010; Thompson et al.,, 2020). We first describe a single species version of the model, before extending this to a multi-species model that also accounts for interactions and competition between coral types.

2.2.1 Generalised logistic growth

We begin with the Richards’ growth model (Richards,, 1959; Tjørve and Tjørve,, 2010; Simpson et al.,, 2022), also known as the generalised logistic growth model, given by

$$\underbrace{\frac{\text{d}C(t)}{\text{d}t}}_{\text{Rate of change in cover}}=\underbrace{\alpha C(t)}_{\text{Coral growth}}\times\underbrace{\frac{1}{\gamma}\left[1-\left(\frac{C(t)}{K}\right)^{\gamma}\right]}_{\text{Competition for available space}},\quad t>t_{0}$$ (1)

where $K\in(0,100]$ (% of area) is the maximum % of area the coral can feasibly cover, i.e., the available space, $C(t)\in(0,100]$ (% of area) is the hard coral cover at time $t$ (years), $\alpha>0$ (1/years) is the intrinsic rate of cover increase and $\gamma>0$ is a non-dimensional generalisation parameter that controls the shape of the sigmoid growth curve. For example, logistic growth is recovered for $\gamma=1$ and Gompertz growth occurs in the limit as $\gamma\to 0$. The effect of $\gamma$ on the recovery curve is shown in Figure 2. Note that Richards’ growth model is usually presented without the $1/\gamma$ factor as it can be absorbed into the definition of the rate parameter, $\alpha$. However, Simpson et al., (2022) show that this leads to non-identifiability in the rate parameter and furthermore renders comparison of $\alpha$ across different values of $\gamma$ to be meaningless.

Figure 2: The effect of the generalisation parameter, $\gamma$, in Richards’ growth model. The logistic growth model corresponds to $\gamma=1$ (blue solid). For $\gamma\to 0$ we obtain the Gompertz growth model with rapid initial recovery and stronger effect of spatial competition. Other values for $\gamma$ result is alternative recovery models. Here, $K=90\%$, $\alpha=0.5\,(\text{Years}^{-1})$, and $c_{0}=5\%$.

2.2.2 Incorporating biphasic recovery

To incorporate biphasic growth into the model, we build on recent developments by Murphy et al., (2022) who consider a piecewise defined population dynamic model,

$$\frac{\text{d}C(t)}{\text{d}t}=\begin{cases}f_{1}(C(t))&\text{if }t\leq T\\ f_{2}(C(t))&\text{if }t>T\end{cases},$$ (2)

where the rate of change of a population, $C(t)$, switches from $f_{1}(C(t))$ to $f_{2}(C(t))$ at the change point $t=T$. It is important to note that the solution to Equation (2) is still continuous at this change point despite the derivative potentially being discontinuous.

To mathematically model the biphasic recovery patterns identified by Warne et al., (2022), we require $f_{1}(C(t))<f_{2}(C(t))$ that is for the time $t\in(t_{0},T]$ the coral recovery rate parameter is lower than when $t\in(T,\infty)$. We achieve this through defining $f_{1}(C(t))=\alpha_{d}f_{2}(C(t))$ where $\alpha_{d}\in(0,1)$ is the recovery scale factor, and $f_{2}(C(t))$ corresponds to the Richards’ growth model (Equation (1)). We then define a duration, $T_{d}>0$, after which the recovery reverts back to Equation (1). Throughout, we refer to this time $T_{d}$ as the relative change point, since it is duration of the first phase relative to the initial time $t_{0}$, that is $T=T_{d}+t_{0}$. Substituting Equation (1) into the general biphasic model of Murphy et al., (2022) (Equation (2)) leads to the following biphasic Richards’ growth model,

$$\frac{\text{d}C(t)}{\text{d}t}=\begin{cases}\dfrac{\alpha_{d}\alpha}{\gamma}C(t)\left[1-\left(\dfrac{C(t)}{K}\right)^{\gamma}\right]&\text{if }t_{0}<t\leq t_{0}+T_{d}\\[10.00002pt] \dfrac{\alpha}{\gamma}C(t)\left[1-\left(\dfrac{C(t)}{K}\right)^{\gamma}\right]&\text{if }t>t_{0}+T_{d}\end{cases}.$$ (3)

Given a known initial condition $C(t_{0})=c_{0}$ and values for the parameters $\boldsymbol{\theta}=[\alpha,\gamma,T_{d},\alpha_{d},K]$, Equation (3) admits the analytical solution,

$$C(t)=\begin{cases}K\left\{1+[(K/C_{0})^{\gamma}-1]\exp{\left(-\alpha_{d}\alpha(t-t_{0})\right)}\right\}^{-1/\gamma}&\text{if }t_{0}<t\leq t_{0}+T_{d}\\ K\left\{1+[(K/C_{d})^{\gamma}-1]\exp{\left(-\alpha(t-T_{d})\right)}\right\}^{-1/\gamma}&\text{if }t>t_{0}+T_{d}\end{cases},$$ (4)

where $C_{d}=C(t_{0}+T_{d})=K\left\{1+[(K/C_{0})^{\gamma}-1]\exp{\left(-\alpha_{d}\alpha T_{d}\right)}\right\}^{-1/\gamma}$.

2.2.3 Accounting for multiple coral family groups

Due to the substantial differences in the recovery rate of the Acroporidae family of corals, it is common to model their recovery distinctly from other hard coral families (Thompson and Dolman,, 2010; Thompson et al.,, 2020; Warne et al.,, 2022). We consider the joint recovery of $M$ coral family groupings with $C_{m}(t)$ denoting the coral cover of family group $m\in[1,2,\ldots,M]$. This leads to a coupled system of $M$ biphasic Richards’ growth model (Equation (3)) with an shared maximum coral cover term to account for competition for space between species,

$$\frac{\text{d}C_{m}(t)}{\text{d}t}=\begin{cases}\dfrac{\alpha_{d,m}\alpha_{m}}{\gamma_{m}}C_{m}(t)\left[1-\left(\dfrac{\sum_{k=1}^{M}C_{k}(t)}{K}\right)^{\gamma_{m}}\right]&\text{if }t_{0}<t\leq t_{0}+T_{d,m}\\[10.00002pt] \dfrac{\alpha_{m}}{\gamma_{m}}C_{m}(t)\left[1-\left(\dfrac{\sum_{k=1}^{M}C_{k}(t)}{K}\right)^{\gamma_{m}}\right]&\text{if }t>t_{0}+T_{d,m}\end{cases},$$ (5)

where $m=1,2,\ldots,M$ and $\boldsymbol{\theta}=[\alpha_{1},\gamma_{1},T_{d,1},\alpha_{d,1},\alpha_{2},\gamma_{2},T_{d,2},\alpha_{d,2},\ldots,\alpha_{M},\gamma_{M},T_{d,M},\alpha_{d,M},K]$. For this work, we only consider the case of $M=1$ or $M=2$. When $M=1$, we consider only the total hard coral cover, and use the analytical solution given in Equation (4). When $M=2$, we use $C_{1}(t)$ to denote the coral cover of the Acroporidae family and $C_{2}(t)$ to denote the coral cover of all other hard corals. See Section 4, Thompson and Dolman, (2010), and Thompson et al., (2020) for discussion around this choice of splitting.

The analytic solution (Equation (4)) for the single species model (Equation (3)) does not extend to the multispecies case (Equation (5)). Therefore, we rely on numerical schemes to solve the coupled system. We apply the Runge-Kutta-Fehlberg method with adaptive time-stepping (Fehlberg,, 1969) as implemented in the R package, deSolve (Soetaert et al.,, 2010).

2.2.4 The observation process model

Let $\mathbf{C}(t)=[C_{1},C_{2},\ldots,C_{M}(t)]^{\text{T}}$ be the vector of coral cover for each of the $M$ modelled family groups. We then model observation noise using a Gaussian distribution, that is the observed cover, $\mathbf{X}(t)$, at time $t$ is given by

$$\mathbf{X}(t)\sim\mathcal{N}(\mathbf{C}(t),\boldsymbol{\Sigma}(t)),$$ (6)

where $\boldsymbol{\Sigma}(t)=\operatorname*{diag}{\{\hat{s}_{m}(t)^{2}/5\}}$ is a diagonal matrix of standard error estimates where $\hat{s}_{m}(t)^{2}$ is the empirical variance in coral cover for family group $m$ at time $t$ taken across the five transects (See Section 2.1). Alternative observation processes could be considered, such as binomial noise, however, the observed heteroscedasticity in the cover data tends to be different than expected for multinomial noise models, that is, the data demonstrates both overdispersion and underdispersion. Using a Gaussian noise approximation allows for a more flexible treatment of the variability that is also computationally efficient.

2.3 Bayesian uncertainty quantification

Inference is performed within a Bayesian framework, that is, for each recovery trajectory $i\in[1,2,\ldots,N]$ we require samples from the joint posterior distribution (Gelman et al.,, 2013) with density,

$$p(\boldsymbol{\theta}^{(i)}\mid\mathcal{D}^{(i)})=\frac{\mathcal{L}(\boldsymbol{\theta}^{(i)};\mathcal{D}^{(i)})p(\boldsymbol{\theta}^{(i)})}{p(\mathcal{D}^{(i)})},$$

where $\mathcal{L}(\boldsymbol{\theta}^{(i)};\mathcal{D}^{(i)})$ is the likelihood of the observed benthic survey data, $\mathcal{D}^{(i)}$, under the model given parameters, $\boldsymbol{\theta}^{(i)}$, $p(\boldsymbol{\theta}^{(i)})$ is the prior probability density and $p(\mathcal{D}^{(i)})$ is the evidence, that is, the marginal likelihood given by

$$p(\mathcal{D}^{(i)})=\int\mathcal{L}(\boldsymbol{\theta}^{(i)};\mathcal{D}^{(i)})p(\boldsymbol{\theta}^{(i)})\textrm{d}\boldsymbol{\theta}^{(i)}.$$

Conditional on the model solution $\mathbf{C}(t_{j};\boldsymbol{\theta}^{(i)})$ (which will be numerical for $M>1$) we obtain the log-likelihood as a multivariate Gaussian around the model,

$$\log\mathcal{L}(\boldsymbol{\theta}^{(i)};\mathcal{D}^{(i)})=-\frac{1}{2}\left\{n^{(i)}\log(2\pi)+\sum_{j=0}^{n^{(i)}}\left[\log{|\boldsymbol{\Sigma}^{(i)}(t_{j})|}+\Delta\mathbf{C}(t_{j};\boldsymbol{\theta}^{(i)})^{\text{T}}\boldsymbol{\Sigma}^{(i)}(t_{j})^{-1}\Delta\mathbf{C}(t_{j};\boldsymbol{\theta}^{(i)})\right]\right\},$$ (7)

where $\Delta\mathbf{C}(t_{j};\boldsymbol{\theta}^{(i)})=\mathbf{X}^{(i)}(t_{j})-\mathbf{C}(t_{j};\boldsymbol{\theta}^{(i)})$ is the difference between the observations and modelled mean recovery at time $t_{j}$ with $j=0,1,\ldots,n^{(i)}$.

We assume independent priors, that is,

$$p(\boldsymbol{\theta}^{(i)})=\prod_{m=1}^{M}p(\alpha_{m}^{(i)})p(\gamma_{m}^{(i)})p(T_{d,m}^{(i)})p(\alpha_{d,m}^{(i)}).$$ (8)

Note that the maximum coral cover, $K^{(i)}$, is not inferred since it is standard practice to estimate this directly from the LMTP and MMP Data (Thompson et al.,, 2020) (See Section 4 for discussion). For the remaining parameters we have, $\alpha_{m}^{(i)}\sim\mathcal{U}(0,1)$, $\alpha_{d,m}^{(i)}\sim\mathcal{U}(0,0.9)$, $\gamma_{m}^{(i)}\sim\mathcal{U}(0,2)$, $T_{d,m}^{(i)}\sim\mathcal{U}(0,D^{(i)})$ where $D^{(i)}=t_{n}^{(i)}-t_{0}^{(i)}$ is the duration of the $i$th recovery trajectory and $\mathcal{U}(a,b)$ denotes a uniform distribution over the interval $[a,b]$. By using uniform priors, we aim for the priors to be vague or uninformative within the physically feasible ranges, however, we note that general guarantees of uninformativity are not practically available (Efron,, 1986; Warne et al.,, 2019).

For our model, the evidence term is intractable. As a result, we can only evaluate the posterior density up to a normalising constant using Equation (7) and Equation (8). To deal with this, we apply Markov chain Monte Carlo (MCMC). In our case, the commonly applied sampler available in the Stan environment (Carpenter et al.,, 2017) or other samplers based on Hamiltonian Monte Carlo (Duane et al.,, 1987) are not appropriate due to a discontinuity that arises in $\nabla\log\mathcal{L}(\boldsymbol{\theta}^{(i)};\mathcal{D}^{(i)})$ from the change-points $T_{d,m}$. While numerical smoothing of the discontinuities is feasible, we avoid this by selecting a more general MCMC sampler that does not require $\log\mathcal{L}(\boldsymbol{\theta}^{(i)};\mathcal{D}^{(i)})$ to be differentiable. Specifically, we use the robust adaptive random-walk Metropolis-Hastings MCMC scheme of Vihola, (2012) that is implemented within the R package, adaptMCMC (Scheidegger,, 2021).

For each trajectory, $\mathcal{D}^{(1)},\mathcal{D}^{(2)},\ldots,\mathcal{D}^{(N)}$, we use four independent chains to assess standard convergence diagnostics of effective sample size (ESS) and the Gelman $\hat{R}$ statistic (Gelman and Rubin,, 1992; Vehtari et al.,, 2021) using the coda package (Plummer et al.,, 2006). Sampling is terminated when $\hat{R}\leq 1.1$ and $\text{ESS}\geq 200$ for each parameter. This results in $N=737$ sets of posterior samples. Thinning is applied for chains with very large ESS to save memory requirements while still ensuring the convergence criteria are satisfied.

To assess the model fit, we sample from the posterior predictive distribution with density,

$$p(\mathcal{D}_{s}^{(i)}\mid\mathcal{D}^{(i)})=\int p(\mathcal{D}_{s}^{(i)}\mid\boldsymbol{\theta}^{(i)})p(\boldsymbol{\theta}^{(i)}\mid\mathcal{D}^{(i)})\,\text{d}\boldsymbol{\theta}^{(i)},$$ (9)

where $\mathcal{D}_{s}^{(i)}$ is simulated trajectory data using the biphasic recovery model (Equation (5)) and observation process (Equation (6)) given a random draw from the posterior samples as obtained via MCMC. Using these simulations, we obtain a probability density estimate for the model uncertainty accounting for all of the underlying parameter uncertainty. This density is then used to estimate 50%, 90%, 95%, and 99% credible intervals (CrI). This enables the model fit to be assessed by inspecting the location of the observed data within these reported intervals.

3 Results

We calibrate the biphasic generalised recovery model using Bayesian analysis for the $N=737$ recovery trajectories as described Section 2. The resulting output is posterior distributions and posterior predictive samples for each of the $N=737$ recovery trajectories. In this section, we provide summaries of the overall model fitness in terms of convergence diagnostics and posterior predictive probabilities. We also present example outputs using four representative recovery trajectories (Figure 3): site 2 of Gannett Cay reef between 2001–2009 (Figure 3(a)), site 1 of Lady Musgrave reef between 1992–2003 (Figure 3(b)), and site 1 of Thetford reef between 1994–1999 (Figure 3(c)) and between 2001–2009 (Figure 3(d)). For ease of description we refer to these recovery trajectories by their reef name and start year, that is Gannett Cay 2001 (Figure 3(a)), Lady Musgrave 1992 (Figure 3(b)), Thetford 1994 (Figure 3(c)), and Thetford 2001 (Figure 3(d)).

Figure 3: Example recovery trajectory data used to highlight the efficacy of our model in capturing different recovery patterns observed for hard coral (solid black lines) broken down into the Acroporidae (dotted black lines) and all other hard corals (dashed black lines). Error bars indicate standard errors for the coral cover estimates.

These example recovery trajectories are selected to demonstrate the range of bisphasic recovery patterns that are observed across the GBR while also providing interesting interplay between the single species and two species models.

3.1 MCMC convergence diagnostics

For each recovery trajectory we obtain samples from the joint posterior distribution under both the single-species (Equation (3)) and the two-species (Equation (5)) versions of the biphasic Richards’ growth model. For the single species model we infer the intrinsic growth rate $\alpha$, the shape parameter, $\gamma$, the relative change-point, $T_{d}$, and the scale factor, $\alpha_{d}$. For the two-species model we infer equivalent parameters for the Acroporidae family, $\alpha_{A},\gamma_{A},T_{d,A},\alpha_{d,A}$, for all other hard corals, $\alpha_{C},\gamma_{C},T_{d,C},\alpha_{d,C}$. In both cases, we used standard methods to estimate maximum coral cover for reefs based on abiotic cover and transient silt records (Jonker et al.,, 2008; Thompson and Dolman,, 2010; Thompson et al.,, 2020).

As described in Section 2.3 we terminate sampling when four independent chains achieve $\widehat{R}\leq 1.1$ and ESS $\geq 200$. Of the $N=737$ recovery trajectories, $38$ did not meet this convergence criteria after 200,000 MCMC iterations. Most of these failed convergences where due to ESS $<200$, with only 11 failing due to $\hat{R}>1.1$. The lower ESS indicates slower mixing in the MCMC chains and a likely cause is the parameter non-identifiability issues that arises when there a very weak biphasic recovery pattern leading to a trade-off between $T_{d}$, $\alpha_{d}$ and $\gamma$, that is the data does not have enough information to distinguish between a single-phase model with $T_{d}=0$ and a biphasic model $T_{d}>0$. In the very few cases where $\hat{R}>1.1$, it seems that the proposed model is a poor fit, and we discuss potential improvements to account for these cases in the discussion (Section 4). That is, we are able to reliably characterise model parameters for $95\%$ of the available recovery trajectories ($N=699$) with only $1.5\%$ ($N=11$) failing due to serious convergence issues arising from poor model fit.

Figure 4: Scatter plots summarising the MCMC convergence diagnostics for each recovery trajectory after a maximum of 200,000 iterations per chain and four independent chains. Convergence thresholds are indicated: $\hat{R}\leq 1.1$ (vertical red dashed line), $\text{ESS}\geq 200$ (horizontal red dashed line). (a) For the majority of recovery trajectories, the MCMC diagnostics comfortably meet the thresholds. (b) Of the small number of recovery trajectories that fail in their MCMC convergence diagnostics, most are close to the threshold.

Figure 5: Trace plots showing good mixing and agreement between four independent MCMC chains (green, orange, purple, and magenta solid lines) in all parameters. (a) Gannett Cay 2001, (b) Lady Musgrave 1992, (c) Thetford 1994, and (d) Thetford 2001.

Table 1: Convergence diagnostic results for the four example recovery trajectories.

Recovery Trajectory	Parameter	$\hat{R}$ [Upper $95\%$CI]	ESS
	$\alpha$	$1.02$ $[1.06]$	429
Gannett Cay site 2 2001	$\alpha_{d}$	$1.01$ $[1.02]$	1092
	$\gamma$	$1.01$ $[1.04]$	990
	$T_{d}$	$1.02$ $[1.05]$	227
	$\alpha$	$1.00$ $[1.00]$	1074
Lady Musgrave site 1 1992	$\alpha_{d}$	$1.00$ $[1.01]$	707
	$\gamma$	$1.00$ $[1.01]$	446
	$T_{d}$	$1.00$ $[1.01]$	454
	$\alpha$	$1.00$ $[1.01]$	969
Thetford site 1 1994	$\alpha_{d}$	$1.08$ $[1.15]$	287
	$\gamma$	$1.00$ $[1.02]$	1549
	$T_{d}$	$1.01$ $[1.03]$	589
	$\alpha$	$1.01$ $[1.01]$	576
Thetford site 1 2001	$\alpha_{d}$	$1.01$ $[1.03]$	395
	$\gamma$	$1.01$ $[1.04]$	365
	$T_{d}$	$1.01$ $[1.03]$	402

In the context of the four examples trajectories (Figure 3), the convergence diagnostics in the form of MCMC trace plots, parameter-wise ESS and $\widehat{R}$ statistics are shown in Figure 5 and Table 1 for the single species model. The trace plots show good mixing and agreement across all four independent MCMC chains (Figure 5). In addition the estimated $\widehat{R}$ statistics including the upper $95\%$CI are well below the standard threshold of $1.1$ (with exception of the $\widehat{R}$ upper $95\%$CI for $\alpha_{d}$ at Thetford 1994), and the ESS for all parameters exceed 200. Given these results, it is reasonable to conclude our convergence criteria are appropriate and that the Markov chains have reached stationarity.

3.2 Assessment of model fitness

For each of the $N=699$ recovery trajectories in which the MCMC sampling satisfied convergence criteria (Section 3.1), we perform posterior predictive checks. That is, we evaluate the location of the observed coral cover within the posterior predictive distribution. For the model to be representing the data accurately we would expect, for any $0<\beta<100$, that at least $\beta\%$ of the observations are within the $\beta\%$CrI of the relevant posterior predictive distribution.

For each observation, we calculate the smallest $\beta$ such that the observation is within $\beta\%$CrI. To do this, samples are generated for the posterior predictive distribution (Equation (9)) and the observed quantiles for each coral cover observation is calculated using,

$$Q_{\textrm{obs}}^{i,j}=\mathbb{P}(\mathbf{X}_{s}^{(i)}(t_{j})\leq\mathbf{X}^{(i)}(t_{j})\mid\mathcal{D}^{(i)})=\int_{0}^{\mathbf{X}^{(i)}(t_{j})}p(\mathbf{X}_{s}^{(i)}(t_{j})\mid\mathcal{D}^{(i)})\,\textrm{d}\mathbf{X}_{s}^{(i)}(t_{j}),$$ (10)

where $\mathbf{X}_{s}^{(i)}(t_{j})$ and $\mathbf{X}^{(i)}(t_{j})$ are, respectively, the simulated and observed coral cover for the $i$th recovery trajectory at the $j$th time point. Given the observed quantile, then the smallest $\beta$ for the CrI containing the observation is given by $\beta^{i,j}=100(1-2Q_{\textrm{obs}}^{i,j})$ when $Q_{\textrm{obs}}^{i,j}<0.5$, and $\beta^{i,j}=100(2Q_{\textrm{obs}}^{i,j}-1)$ otherwise. Hence, we can estimate the proportion of observations that within the posterior predictive $\beta\%$CrI using

$$\hat{p}_{\beta}=\frac{\sum_{i=1}^{N}\sum_{j=1}^{n^{(i)}}\mathbb{I}\left(\beta^{i,j}<\beta\right)}{N\sum_{i=1}^{N}n^{(i)}}.$$ (11)

Here, $\mathbb{I}\left(\beta^{i,j}<\beta\right)=1$ if $\beta^{i,j}<\beta$ and $\mathbb{I}\left(\beta^{i,j}<\beta\right)=0$ otherwise. The standard error for the estimate of is given by

$$\hat{s}_{\beta}=\sqrt{\frac{\hat{p}_{\beta}(1-\hat{p}_{\beta})}{N\sum_{i=1}^{N}n^{(i)}}}.$$ (12)

Comparing estimates of $\hat{p}_{\beta}$ with $\beta=1,2,\ldots 99$ provides a summary of how well posterior predictive distributions fit the data. As demonstrated in Figure 6(a), we observe that for $\beta\in\{1,91\}$, we have $\beta<\hat{p}_{\beta}\pm 1.96\hat{s}_{\beta}$, indicating the data is well represented by the bulk of the posterior predictive distributions. When $\beta\in\left\{92,93,\ldots,97\right\}$ we observe $\hat{p}_{\beta}-1.96\hat{s}_{\beta}<\beta<\hat{p}_{\beta}+1.96\hat{s}_{\beta}$, indicating the posterior predictive distributions are also tracking the more extreme data points at the expected frequencies. In the extreme tails, $\beta\in\left\{98,99\right\}$ we have $\beta>\hat{p}_{\beta}\pm 1.96\hat{s}_{\beta}$ and conclude the very extreme variations in the data are under-represented, through the margin of inaccuracy here is still only a few percent. We conclude that overall the performance of the biphasic Richards’ model in capturing the variation and dynamics of coral reef recovery is extremely high.

Figure 6: Posterior predictive analysis over the entire GBR data set including $N=699$ recovery trajectories for which the MCMC sampling converged according to our criteria. (a) The estimated proportion of data that falls within the $\beta\%$CrI, $\hat{p}_{\beta}$, (solid black lines, with error bars for $95\%$ confidence intervals) for each level of $\beta$ (red line for $\beta=\hat{p}_{\beta}$). (b) Bivariate density plot of the smallest containing posterior predictive CrI for each coral cover observation (yellow indicates high density and dark blue low density).

In Figure 6(b) a density plot shows the distribution of smallest containing CrIs for each coral cover observation. For coral covers below $25\%$ cover the pattern, the bulk of smallest containing CrIs are well below the $25\%$CrI and the distribution skewed to the right, indicating high posterior probability density for most of the observed data. As coral cover increases, especially above $40\%$ cover, the distribution of smallest containing CrIs is more diffuse. From this we conclude that the majority of the data overdispersion is occurring in regions of higher cover. Importantly, the greatest accuracy in the model is for cover below $30\%$, which is where the majority of the GBR sits. It is also important to emphasise that at higher cover we still observe good fitness, but the data is more variable in the extremes than the model predicts.

Figure 7(a)–(d) demonstrates example posterior predictive distributions for both the single species model (Equation (3)), and the two species model (Equation (5)) that describes the competition for space between fast and slow growing corals. The posterior predictive plots are compared with the observed data in both the single-species and two-species representations for each of the example recovery trajectories: Gannett Cay 2001 (Figure 7(a)), Lady Musgrave 1992 (Figure 7(b)), Thetford 1994 (Figure 7(c)), and Thetford 2001 (Figure 7(d)). In each case, there is excellent agreement between the model predictions and the observed data with majority of data points lying within the 95%CrI of the posterior predictive distribution.

Figure 7: Posterior predictive distributions for example recovery trajectories. (a)–(d) The observed data is given for the total cover for hard coral (solid black lines) broken down into Acroporidae (dotted black lines), and other hard corals (dashed black lines); and model predictions with credible intervals (50%, 90%, 95% and 99% CrI) as shaded regions for total hard coral (green), Acroporidae (red), and other hard corals (blue).

3.3 Parameter inference

For a given recovery trajectory, marginal posterior densities for the model parameters are estimated directly from the MCMC samples obtained by targeting the full posterior distribution (Sections 2.3 and 3.1). These marginal densities can be used to obtain insight into the underlying mechanisms leading to a particular recovery pattern for a recovery trajectory of interest. In Figure 8, univariate posterior marginal probability densities are shown by coral family group for each of the examples recovery trajectories: Gannett Cay 2001 (Figure 8(a),(e),(i),(m)), Lady Musgrave 1992  (Figure 8(b),(f),(j),(n)), Thetford 1994 (Figure 8(c),(g),(k),(o)), and Thetford 1994 (Figure 8(d),(h),(l),(p)). Among these examples, our model identifies a range of different recovery patterns that could be of interest to reef managers.

Across each of the example trajectories there are some common features in the marginal posterior distributions for the recovery rates (Figure 8(a)-(d)). Firstly, the two species model tends to indicate the recovery rate for the Acroporidae family is greater than the other hard corals, and the recovery rate for the single species model tends to sit between the two sub-groups with the relative weighting between the two subgroups indicative of their relative coral cover contribution (Figure 7 and Figure 8(a)-(d)). This is expected since members of the Acroporidae family are known to be fast growing by comparison to most other hard corals (Thompson and Dolman,, 2010; Ortiz et al.,, 2021). An interesting deviation from this pattern is in the estimates for the recovery rate for Thetford 2001 (Figure 8(d)) where the total hard coral estimates are more dispersed than either of the sub-groups. An explanation for this could be the more complex dynamics in the other hard corals with an initial decline before an increase in cover coinciding with the Acroporidae corals exceeding the cover of the other hard corals (Figure 7(d)). This results in a more complex recovery pattern for the total hard coral cover.

The shape parameter in the Richards’ growth model shows a variety of patterns across the example recovery trajectories (Figure 8(e)–(h)) with some examples approximating Gompertz growth ($\gamma\to 0$), others approximating logistic growth ($\gamma=1$), and others tending toward other forms from Richards’ growth ($\gamma>1$). This is not surprising given the diversity in coral species in the GBR. It is worth noting that the posterior distributions for the shape parameter in both Thetford 1994 (Figure 8(g)) and Thetford 2001 (Figure 8(h)) are almost identical for the total and other hard corals, which would be expected for the same reef. The difference observed in the posterior for the shape parameter for the Acroporidae is likely due to the data containing little information on the sigmoid shape since the recovery for Acroporidae is near linear (Figure 7(d)). An alternate explanation could be a change in the specific community composition in the Acroporidae group at Thetford, we do not explore this further here. We highlight, however, that identifying similarities and differences in the recovery patterns of a specific reef over time is a key use case for coral recovery models in marine ecology and monitoring (Thompson et al.,, 2020; Ortiz et al.,, 2018).

Figure 8: Univariate marginal posterior densities for (a)–(d) the recovery rate $\alpha$, (e)–(h) shape parameter $\gamma$, (i)–(l) change-time $T_{d}$, and (m)–(p) scale factor $\alpha_{d}$. Densities are given for total hard coral (green) broken down into Acroporidae (red), and other hard corals (blue). Results are shown for the four examples recovery trajectories: (a),(e),(i),(m) Gannett Cay 2001, (b),(f),(j),(n) Lady Musgrave 1992, (c),(g),(k),(o) Thetford 1994, and (d),(h),(l),(p) Thetford 2001.

For the parameters that relate to the biphasic recovery, the example recovery trajectories demonstrate the range of different behaviours between different locations and times, and between the single species and two species representations. In the following paragraphs we highlight key features in each example. One important property of our biphasic model (Equation (3) and Equation (5)) is that as the scale factor $\alpha_{d}\to 1$, the change-time becomes structurally non-identifiable. Thus whenever the posterior density for the change-time is very flat, and the posterior density for the scale factor skewed to the left with bulk above $0.6$, then we conclude there is no evidence for biphasic recovery.

In the analysis of Gannett Cay 2001, the single species model indicates little evidence of a biphasic recovery pattern for total hard coral cover. This conclusion is drawn from the relatively flat posterior for the change time in the (Figure 8(i)), and the bulk of the scale factor is tending toward larger values (Figure 8(m)). However, when the data is interpreted with the two species model, there is a strong signal for biphasic recovery patterns are present in the Acroporidae group with strong posterior support for change-time between $1$ to $4$ years (Figure 8(i)) and a small value for the scale factor (Figure 8(m)). For the other hard corals, the posteriors for the change-time and scale factor are very similar to the single species model with little evidence for biphasic recovery.

At Lady Musgrave 1992, there is a strong indication for biphasic recovery patterns in total hard coral cover with a change-time of between $1$ to $3$ years (Figure 8(j)) and scale factor less than $0.6$ (Figure 8(n)). When using the two species model, it is clear that a biphasic recovery pattern in the Acroporidae family, with a change-time of around $1.5$ to $2.5$ years and scale factor less that $0.5$, provides the main driver for the overall biphasic recovery for the total hard coral cover. This is expected at Lady Musgrave where the Acroporidae family accounts for almost 100% of corals (Figure 7(b)) (Simpson et al.,, 2022; Murphy et al.,, 2022; Simpson et al.,, 2023). However, the uncertainty in the biphasic parameters for the small population of other hard corals in the two species model (Figure 7(b)) is reflected in the larger variance in the change-time and scale factor posterior densities for the total hard corals (Figure 8(j),(n)).

Analysis for Thetford 1994 and 2001, are consistent with each other in relation to biphasic recovery. In both recovery trajectories, there is evidence for biphasic recovery patterns in the total hard coral cover (Figure 8(k)–(l),(o)–(p)) (though there is substantially more uncertainty in the results following 1994). Just as in the analysis of Lady Musgrave 1992, one of the species in the two species model seems to be the main driver for biphasic recovery pattern in the total hard coral cover, however, for Thetford it is the other hard corals rather than the Acroporidae. We arrive at this conclusion by noting, the change-time and scale factor posteriors densities for the other hard corals follow qualitatively a similar pattern to the equivalent posteriors for total hard corals (Figure 8(k)–(l),(o)–(p)). The Thetford 2001 recovery trajectory is more informative in relation to the biphasic recovery pattern (Figure 8(k),(o)), however, both Thetford 1994 and 2001 are qualitatively similar with no evidence for biphasic recovery patterns in the Acroporidae family and agreement between other hard corals and total hard corals in support of biphasic recovery patterns (Figure 8(k)–(l),(o)–(p)). Due to the consistency between distinct trajectories, this would support the theory that this biphasic recovery pattern is a feature of the environment or ecosystem at Thetford reef.

4 Discussion

In this work, we perform Bayesian parameter inference and posterior predictive checks for $N=699$ recovery trajectories across the GBR. We find both the single species and two species biphasic Richards’ growth model is flexible enough to obtain good model posterior predictive performance. This is impressive in of itself given the range of possible recovery patterns that could arise from the complex interaction of different benthic species. Using four recovery trajectories as exemplars, we highlight the utility of the method for uncertainty quantification in key model parameters related to the Richards’ growth shape, recovery rate, and biphasic recovery patterns. Interpretation of these parameters both from a single species and two species viewpoint could support reef management and interventions following significant disturbance events.

Accurately predicting and interpreting recovery patterns for diverse locations across the GBR is critical for management to support reef recovery, resistance and adaptation into the future (Ortiz et al.,, 2021; Thompson et al.,, 2020). Motivated by prior work demonstrating potential widespread deviations from classical Gompertz or logistic recovery models (Warne et al.,, 2022; Murphy et al.,, 2022; Ortiz et al.,, 2018), we develop a biphasic extension to the multispecies Richards’ growth model. This framework is flexible enough to capture a range of possible sigmoid recovery patterns and biphasic recovery patterns where a period of inhibited recovery occurs in the year immediately following the disturbance event (Warne et al.,, 2022). In this paper, we formulate this new model and demonstrate its efficacy in accurately capturing the recovery patterns of nearly every recovery trajectory available within LTMP reef monitoring data, which is spatially and temporally the largest reef monitoring dataset in the world. Therefore, our modelling approach is describing well the recovery patterns observed in the largest and most interconnected coral reef worldwide. Through Bayesian uncertainty quantification we can reliably characterise the combination of parameters that lead to particular recovery behaviours for a reef and provide forecasts of coral recovery following a disturbance event.

In this work, we considered the LTMP data from a single species perspective, based on total hard coral cover, and from a two species perspective that separates the fast-growing Acroporidae coral families from the other hard corals (Thompson and Dolman,, 2010; Thompson et al.,, 2020). We believe that both models are important for management. Firstly, the single species model enables a macroscopic view of the reef recovery patterns that is useful as a reef health indicator (Thompson et al.,, 2020) and the notion of a biphasic recovery has a direct interpretation in terms of expected future coral cover that, in turn, can support intervention triage (Warne et al.,, 2022). In addition, the single species model is extremely efficient for the purposes of model calibration due to the lower dimensionality of the parameter space and the analytic solution to the biphasic Richards’ growth model (Tjørve and Tjørve,, 2010; Simpson et al.,, 2022). However, the two species model provides more insight around the potential mechanisms driving any biphasic recoveries that are observed using the single species model. This provides additional information to management authorities regarding the types of interventions that may be applicable. Unfortunately, this additional granularity of interpretation comes with a computational cost as numerical schemes are needed to solve the governing equations and there are more parameters that can lead to slower convergence in the MCMC sampling. Therefore, we recommend a two-stage process for monitoring in which the single species model is used to triage reefs affected by biphasic recovery, then the two species model is introduce to assist in planning interventions. Future work could also explore the use of frequentist methods based on likelihood profiles (Simpson et al.,, 2023) as a preconditioning step to accelerate the MCMC sampling and reduce the computational burden as proposed by Warne et al., (2024). Such improvements would be important if further refinements are considered, such as the inclusion of additional hard coral family divisions or the inclusion of other benthic life forms that compete for space with coral such as macro-algae and sponges (Thompson and Dolman,, 2010; Mumby,, 2009), or the inclusion of spatial hierarchical patterns (Johns et al.,, 2014).

One limitation to our process is the pre-estimation of maximum coral cover. While this is a standard procedure in the coral monitoring literature (Thompson et al.,, 2016, 2020), there is evidence to suggest that these estimates could be overestimates to the effective maximum coral cover (Cresswell et al.,, 2024). In addition, it is also mathematically problematic to pre-estimate long-term phenomena, in general. We also note that this pre-estimation of the maximum coral cover is one of the main reasons for the relatively few situations in which we obtain poor model fits or convergence issues with MCMC sampling. Therefore, future work should also consider the inclusion maximum coral cover parameter in the model calibration process. This will be especially important if there is latent competition with other species, like macro-algae, that are not currently explicitly modelled.

We have only considered deterministic models designed to capture the average recovery patterns. As a result, the noise is assumed to be entirely due to the monitoring process (Jonker et al.,, 2008; Thompson et al.,, 2016). However, at a colony level, reef growth and competition with other species is a stochastic process (Álvarez–Noriega et al.,, 2023; Bozec et al.,, 2021; McDonald et al.,, 2023). The net effect of this intrinsic noise could be captured at the reef scale using stochastic differential equation (SDE) models. While the use of SDEs generally requires more complex computational approaches to perform parameter inference (Cranmer et al.,, 2020; Drovandi and Pettitt,, 2010; Sisson et al.,, 2018; Wang et al.,, 2024; Warne et al.,, 2020, 2024), there are known benefits in using SDE models when intrinsic noise is relevant, as they lead to more accurate parameter inferences and can resolve issue related to parameter non-identifiability (Browning et al.,, 2020; Munsky et al.,, 2009). This could be of considerable benefit, as practical identifiability issues arise in the uncertainty quantification analysis for some parameters (Figure 8) and is it not feasible to collect more data.

Prior work by Warne et al., (2022) demonstrated the importance of taking biphasic recovery patterns in coral reefs into account. Through the application of our Bayesian framework that accounts for this phenomena, the construction of coral reef health indicators that are robust to deviations in standard assumptions could be constructed (Simpson et al.,, 2022; Thompson and Dolman,, 2010; Thompson et al.,, 2020). This new approach enables a more nuanced view of expected reef recovery and provides a modelling basis to investigate the biological, ecological and environmental factors that contribute to biphasic recovery patterns (Warne et al.,, 2022). New measures of reef decline can also be established with this approach. For example, tracking the distribution of the change time and scale factor for a given reef over multiple recoveries can identify steady increases in the duration of the slower initial recovery phase that could indicate degradation in reef health. We anticipate our modelling and analysis approach to have direct benefit for the management of the GBR. Finally, through the use of decision-theoretic approaches to sampling design (Akinlotan et al.,, 2024), understanding biphasic recovery patterns could enable new monitoring protocols that maximise information for decision makers.

Acknowledgements

This work was supported through the Australian Institute of Marine Science (AIMS) and the Centre for Data Science at the Queensland University of Technology. Computational resources were provided by the eResearch Office, Queensland University of Technology.