Following a disturbance on a coral reef, and in the absence of further extreme environmental conditions, the change in average precentage coral cover (P) through time (Fig. 1) often can be approximated by the logistic growth equation. To add environmental complexity and stochastic extremes we considered that the physiological response of corals depends on temperature stress, which in turn depends on irradiance32. A newly constructed hybrid stochastic model, given below, governs the dynamics of that response and the occurrences of extreme-temperature events:

$$\frac{dP}{dt}=rP(1-\frac{P}{K})-\gamma TP\,{\rm{for}}\,t

otin ({t}_{k},\,{t}_{k}+h)\,{\rm{and}}\,{t}_{0}\le t\le {t}_{N}{\rm{with}}\,P({t}_{0})={P}_{0},$$ (1a)

$$P(t)=P({t}_{k}){e}^{f(\varepsilon )t},\,{t}_{k}\le t\le {t}_{k}+h\,{\rm{for}}\,k=1,2,\ldots ,\,n\le N,$$ (1b)

where r is the intrinsic rate of increase in coral cover per year, K is the steady-state coral cover equilibrium point (%), T is the sea-surface temperature (°C), and where t k ’s are the years when extreme thermal-stress events occur. N is an integer that indicates the number of years we wish to run the simulations, and in our case it was through to the year 2100. In addition, h is the duration of a thermal-stress anomaly (weeks), which was extracted from a truncated gamma distribution with a minimum of 2 weeks and a maximum of 8 weeks, and γ is a coefficient that impacts coral populations through temperature changes, which affects percentage coral cover. During an extreme-temperature event, coral cover is modeled by Equation 1b, where ε, which we call an ‘extreme-event coefficient’, is the sudden increase in temperature, herein governed by short-term climate oscillations such as El Niño events, which result in mass coral bleaching and mortality. The function f(ε) controls the dynamics, specific to each coral species or coral assemblage, depending on sudden temperature increases (measured by ε). Since the coral populations decrease during temperature anomalies, f(ε) is a negative function for all times. Using the present biological data the function was estimated as −ε3. These values would most likely vary geographically and for different coral assemblages, and could vary into the future, with anomalous temperature events not necessarily coupled with El Niño cycles5, although the generic construct of the model will remain the same.

Figure 1 Coral-population data for the Scott Reef system in north Western Australia (Gilmour et al.15), and for Sesoko Island, southern Japan (van Woesik et al.14). (a) Average coral cover (%). (b) Average coral cover (%) of Acropora spp. and Porites spp. from Gilmour et al.15. Gray bars indicate the timing of the El Niño-driven mass-bleaching and mortality event between 1997 and 1998. Note that the y-axes differ between panels. Full size image

Sea-surface temperature (in Celsius), T(t), was obtained by curve fitting satellite data (see below), using the following equation:

$$T(t)={I}_{ave}({a}_{1}\,\cos (t)+{a}_{2}\,\sin (t))+\lambda t+{a}_{3},$$ (2)

where I ave is the annual average of irradiance (in photosynthetic available radiation, PAR), and the parameters a 1 and a 2 are the rates at which irradiance affects the change in temperature (or the rate of change in Celsius per PAR) (in Celsius/PAR). We assume that seasonal irradiance does not vary among years. Here, λ is defined as the ‘climate-change coefficient’ (Celsius time−1), where λt controls the temperature increase over time, in years, and a 3 is the annual average of sea surface temperature (in Celsius) for λ = 0 under normal, non-anomalous years (Table 1). To estimate these parameters we used nonlinear optimization using the program Mathematica®.

Table 1 Estimated coefficients of sea-surface temperature functions in Equation 2 for each Representative Concentration Pathway (RCP) scenario (IPCC 2013). Full size table

Biological data

To estimate the parameters in Equation 1, we first estimated recovery rate (r) and carrying capacity (K). We obtained data sets of coral cover from two geographic regions: Sesoko Island, Okinawa, in southern Japan11,14, and Scott Reef in Western Australia15. Total coral cover estimates were obtained from Sesoko Island from 1997–201011,14. The corals at Sesoko Island experienced an anomalous thermal-stress event in 1998, which was driven by the El Niño Southern Oscillation, which reduced coral populations significantly11. The Sesoko Island reefs also experienced a milder thermal-stress event in 2001 (Fig. 1). In Western Australia, Gilmour et al., see ref.15, collected data on total cover, and the cover of dominant reef-building genera, including Acropora and Porites, from Scott Reef from 1994–2010 (Fig. 1). The reefs also suffered coral losses during the 1998 anomalous thermal-stress event (Fig. 1).

We could not, however, with any degree of certainty, simultaneously fit all the unknown parameters in Equation 1a and also predict coral cover into the future. Therefore, we used a two-step process. First, we estimated r, K and γ as inverse problems using a Bayesian platform. The values from the Bayesian posterior distributions were then fitted to parametric distributions. Second, we used these best-fit distributions to predict the response of coral populations to likely thermal-stress scenarios (see below) through to the year 2100. We estimated r and K values (Equation 1a) at each study locality as an inverse problem using the logistic growth function within a Bayesian framework using a lag-1 temporal auto-regressive error structure on the residuals (after ref.34). The successive r-values for the different localities were derived from a sample size of 50,000 from the posterior distribution, estimated using OpenBUGS. These r-values, from the posterior distribution, were low, <0.4, and were therefore fitted to Beta distributions using the R package ‘fitdistrplus’35. Random samples were taken from these estimated Beta distributions for every time step in the predictive model. We used a similar process to estimate γ-values using Equation 3 below. Although temperatures change through the seasons, in the tropics they only range about 3 °C, which we assumed was small enough to fix T, and thereby use mean temperature and the solution of Equation 1a (i.e., Equation 3) to estimate γ-values. Therefore, to derive a meaningful estimate of γ, we treated temperature,T, and K as constants, and treated r as a stochastic parameter, using the following:

$$P(t)=\frac{(r-\gamma T){P}_{o}}{\frac{r}{K}\,{P}_{o}(1-ex{p}^{-(r-\gamma T)t})+(r-\gamma T)\,ex{p}^{-(r-\gamma T)t}}.$$ (3)

Temperatures, Thermal stress, El Niño, and Climate Change

Sea-surface temperature (SST) data were obtained for each reef from the Moderate-resolution Imaging Spectroradiometer (MODIS) on board the Aqua satellite platform that was launched in 200236. Average monthly MODIS Aqua Global Level 3 mapped mid-infrared nighttime SSTs from January 2003 (the earliest available date) through December 2010 were obtained at a 4.6 km spatial resolution for each locality (Sesoko Island 26.65854°N 127.8611°E; Scott Reef 14.167959°S, 121.87968°E). Monthly SST records were extracted from each study site and mean monthly SST temperature was calculated.

To estimate the distribution of the frequency of thermal-stress events, we obtained Oceanic Niño Index (ONI) data from 1950–2016. The ONI is the standard metric used to identify El Niño (warm) and La Niña (cool) events in the tropical Pacific, and is a running 3-month average sea-surface temperature anomaly for the Niño 3.4 region37 (5°N–5°S and 120°W–170°W), based on the Extended Reconstructed Sea Surface Temperature (ERSST v4) dataset from the National Oceanic and Atmospheric Administration’s National Climate Data Center. These data were derived from the International Comprehensive Ocean-Atmosphere Dataset at a 2° by 2° spatial resolution (https://www.ncdc.noaa.gov). El Niño events were defined as five consecutive three-month periods, at or above the 0.5 °C anomaly. We categorized El Niño events into one of four strengths based on the magnitude of the anomaly: (1) weak events were characterized as 0.5–0.9 °C anomalies, (2) moderate events were defined as 1.0–1.4 °C anomalies, (3) strong events were defined as 1.5–1.9 °C anomalies, and (4) very strong events were greater than or equal to 2.0 °C anomalies (Table 2). The frequency of thermal-stress events was then estimated for each of four categories: (1) all events; (2) only moderate, strong, and very strong events; (3) only strong and very strong events; and (4) only very strong events (Fig. 2).

Table 2 Average and median periodicity of El Niño events from 1950–2016, based on data from the Oceanic Niño Index (https://www.ncdc.noaa.gov). Full size table

Figure 2 History of the frequency of El Niño events from 1950–2016, where ‘All’ indicates all events, ‘MSV’ indicates moderate, strong, and very strong, ‘SV’ indicates strong and very strong, and ‘V’ indicates very strong events. Full size image

To estimate the distribution of the frequency of thermal-stress events of each category of El Niño, the number of years in between the start of each event was calculated. The mean, median, and standard deviation of the period in between the anomalous-thermal events were also calculated (Table 2). The frequency of extreme events in the model simulations was drawn from a Poisson distribution where the rate parameter was estimated from the mean of the frequency of thermal-stress events of categories 1–3 (i.e., between 3 and 9 years). The historical record shows that very strong thermal-stress events occurred every 21.7 ± 9 years (median = 18) (Fig. 2); therefore very strong events were excluded from this study because the infrequent return period lowers the confidence in these events for our relatively short timeframe (to the year 2100).

It is expected that, by 2100, sea-surface temperatures will increase by approximately 1.8 °C, 2.2 °C and 3.7 °C for the Representative Concentration Pathways (RCPs) 4.5, 6 and 8.5 respectively38. Therefore, the maximum temperature anomaly was set to 3.7 °C bove the mean background temperatures in the models38. Based on the RCP predictions, the corresponding climate-change coefficient values, λ, in Equation 2 were re-calculated using curve fitting (Table 1).

Model implementation