Past the Precipice? Projected Coral Habitability Under Global Heating.

Coral reefs are rapidly declining due to local environmental degradation and global climate change. In particular, corals are vulnerable to ocean heating. Anomalously hot sea surface temperatures (SSTs) create conditions for severe bleaching or direct thermal death. We use SST observations and CMIP6 model SST to project thermal conditions at reef locations at a resolution of 1 km, a 16-fold improvement over prior studies, under four climate emissions scenarios. We use a novel statistical downscaling method which is significantly more skillful than the standard method, especially at near-coastal pixels where many reefs are found. For each location we present projections of thermal departure (TD, the date after which a location with steadily increasing heat exceeds a given thermal metric) for severe bleaching recurs every 5 years (TD5Y) and every 10 years (TD10Y), accounting for a range of post-bleaching reef recovery/degradation. As of 2021, we find that over 91% and 79% of 1 km2 reefs have exceeded TD10Y and TD5Y, respectively, suggesting that widespread long-term coral degradation is no longer avoidable. We project 99% of 1 km2 reefs to exceed TD5Y by 2034, 2036, and 2040 under SSP5-8.5, SSP3-7.0, and SSP2-4.5 respectively. We project that 2%-5% of reef locations remain below TD5Y at 1.5°C of mean global heating, but 0% remain at 2.0°C. These results demonstrate the importance of further improving ecological projection capacity for climate-vulnerable marine and terrestrial species and ecosystems, including identifying refugia and guiding conservation efforts. Ultimately, saving coral reefs will require rapidly reducing and eliminating greenhouse gas emissions.

Introduction

Coral reefs are among the most biodiverse ecosystems on the planet (Veron, 1995). However, over the last decade there has been a rapid global decline in coral health and coral cover due to both local environmental degradation (from destructive fishing practices, overfishing, coastal development, sedimentation, nutrient over‐enrichment, and chemical pollutants, and other causes) and global climate change (increasing ocean heat, sea levels, and ocean acidification; De'ath et al., 2012; Hughes et al., 2017).

Although regional bleaching events had been occasionally observed throughout the twentieth century (Yonge, 1930), the first mass event occurred during the 1982–1983 El Niño. It included effects across the Indo‐Pacific (Coffroth et al., 1990) and was likely more widespread than documented. The first global bleaching event occurred during the 1997–1998 El Niño (Hoegh‐Guldberg et al., 2017). The next global event occurred in 2010, and the third began in 2014 and lasted three years. Over recent decades, 33%–50% of coral reefs have been largely or completely degraded (The International Society for Reef Studies, 2015). Overall, there is great concern about the current state of reefs and for their future, as humans continue to heat the planet (Langlais et al., 2017).

Several prior studies have used sea surface temperature (SST) outputs from global Earth system and climate models (hereafter global models) to assess future bleaching risk (Donner, 2009; Frieler et al., 2013; Hoegh‐Guldberg, 1999; Schleussner et al., 2016; Van Hooidonk et al., 2013, 2016). These studies most often report TD5Y, the year after which a thermal threshold is subsequently surpassed at least once per five years, at GM‐like spatial resolution of ∼100 km2. Severe bleaching projections could better inform local conservation decisions if they could capture spatial structure at ∼1 km (Van Hooidonk et al., 2016). Downscaling global model SST projections can therefore better inform decision‐making, and statistical downscaling compares well to more computationally expensive dynamical downscaling (Van Hooidonk et al., 2015). Here, we provide projections of thermal severe bleaching from an ensemble of CMIP6 global models, and at a spatial resolution of 1 km. After submitting our study and while it was undergoing peer review, another study was published independently that also projects coral futures using CMIP6 models downscaled to 1 km resolution (Dixon et al., 2022); compared to this study, ours uses a more advanced downscaling methodology, considers multiple thermal thresholds that explicitly reference the climatological baseline, and provides projections of departure years under multiple emissions scenarios in addition to projections referenced to Global mean surface temperature anomalies (GMSTA).

Data and Methods

CMIP6 Model Data

We included in the analysis one run (or “member”) from every CMIP6 model available as of 25 December 2021 with monthly SST output for the historical experiment and the four future emissions scenarios SSP1‐2.6, SSP2‐4.5, SSP3‐7.0, and SSP5‐8.5 (SSP is “Shared Socioeconomic Pathway,” O’Neill et al. (2014)). These four scenarios span a range of possible collective human futures in terms of greenhouse gas (GHG) emissions, in order of increasing cumulative emissions, with SSP585 being the highest; the final two digits provide the estimated radiative forcing in 2100 in W/m2. In what follows, we omit the punctuation in the emissions scenario labels. In all, the analysis included 35 members from 35 model groups. The model member chosen was the one with the most experiments run, with ties chosen alphabetically (e.g., “r1i1p1f1” over “r2i1p1f1”). We decided to use only one model member per model group in order to avoid multiple members from a single group from potentially biasing the ensemble mean (In the Supporting Information S1 we present results from a different ensemble with 127 members from 27 groups.) The CMIP6 historical experiment begins in January 1870 and runs to December 2014, while the SSP experiments start in January 2014 and run until at least 2100. We regridded all models to be on the same 1° grid and homogenized all time dimensions to the same mid‐month values. The few models that ran beyond December 2099 were truncated to that month.

Global mean surface temperature anomalies (GMSTA) were estimated using 2 m surface temperatures from 33 global models (available as of 28 August 2020), one member from each of 33 model groups, which were each regridded to the same uniform 1° grid. The area‐weighted mean was taken for each model, and then the mean over every model per scenario was taken. GMSTA were calculated relative to an 1880–1900 baseline.

Observational Data

For performing statistical downscaling and for performing degree heating week estimates at 1 km scale, we use NASA/JPL Multiscale Ultrahigh Resolution (MUR) observational SST data from remote sensing, a 0.01°(∼1 km in the domain of our analysis) gridded daily satellite product, available from 2002 to the present, which increases feature resolution over existing SST analysis products with resolutions of 10–100 km. We average the daily MUR product into a monthly product.

The RMS difference between MUR and the quarter‐degree‐gridded Group for High Resolution Sea Surface Temperature Multi‐product Ensemble median SST analysis is 0.36°C in non‐Arctic regions on a daily comparison basis (Chin et al., 2017). Assuming that both SST datasets are unbiased and have equal variance, we can then estimate the error in MUR at one standard deviation to be 0.25°C on a daily basis, or roughly 0.05°C on a monthly basis. This should be thought of as lower bound on the monthly observational SST uncertainty as it excludes potential systematic biases.

To determine the locations of coral reefs in the global ocean, we use a 4 km resolution reef mask from the NOAA Coral Reef Watch thermal history product, v1.0 (Heron et al., 2016), which yields 989,936 1 km reef pixels with the caveat that some 4 km reef pixels may not be fully populated with 1 km reefs. Any 1°coarse pixel that has fewer than 10 global model output values (due e.g., to some models assuming a land pixel and assigning a null value) is excluded from the analysis. This leaves 773,261 1 km reef pixels remaining.

Degree Heating Week Thresholds

DHW is a thermal stress index developed decades ago by Coral Reef Watch (Liu et al., 2003, 2006). At a given location, the maximum monthly mean (MMM) is determined from a climatology (the climatologically hottest month of the year). Then for each day the MMM is subtracted from that day's SST, and if the result is >=1°C (i.e., a degree or more over the MMM) it is accumulated in a 12‐week running sum. According to Coral Reef Watch, significant bleaching in corals is correlated to DHW values >4 DHW, and severe bleaching is likely and significant mortality can be expected above 8 DHW (Coral Reef Watch, n.d.). The original Coral Reef Watch DHW metric requires a 1°C excursion above MMM before it accumulates a daily value into DHW.

Following all of the previous monthly projection studies (see e.g., Van Hooidonk et al. (2016)), we deviate from the Coral Reef Watch definition by not requiring the >=1°C daily excursion above MMM, which cannot be implemented using monthly time series. This allows fair comparisons to those previous monthly projection studies, which we will discuss below. Furthermore, there is evidence that not requiring the >=1°C daily excursion above MMM increases the skill of the DHW metric at predicting bleaching (DeCarlo, 2020; Kim et al., 2019). To calculate an approximate DHW index, we first create a monthly MUR SST climatology from 2003 to 2014, inclusive, which determines a MMM value at each 1 km coral pixel. We subtract this MMM from the SST time series at that pixel, setting any negative values to zero, and multiply by 4.34 to convert from months to weeks. We then calculate a 3 month running sum, producing a monthly time series of DHW estimates. In what follows, we will use “DHW” to also indicate units of °C‐weeks.

The original Coral Reef Watch 8 DHW severe bleaching threshold is based on a climatology comprised of the seven‐year period of 1985–1990 plus 1993 which excludes SST retrievals compromised by the Pinatubo eruption (Heron et al., 2014), the mean of which is 1988.3. In 2015, Coral Reef Watch updated their DHW product, shifting to a new climatological reference period centered at 1998.5 (Liu et al., 2014). However, as mentioned above, the MUR SST climatology central year is 2008.5. In the two decades spanning these three climatological references, SST in coral‐reef‐containing waters increased by 0.25°C due to anthropogenic global heating, as estimated from the mean of all 1‐degree‐resolution HadISST (an observational SST record, Rayner et al. (2003); National Center for Atmospheric Research Staff (Eds) (n.d.)) grid cells containing coral reef locations, with a 10‐year running mean applied to the resulting time series.

The effect of this anthropogenic increase in the climatological baseline is often neglected, but it has a critical impact on DHW metrics. We empirically determined the (linear) relationship between the climatological central year and the DHW threshold required to keep departure year projection estimates constant (see Supporting Information S1 for the detailed methodology). Using subscripts to denote the integer part of the climatological central years discussed above, we found that, for example, in other words, fully specifying a DHW threshold requires two numbers, the threshold and the climatological center year used to calculate it; and an 8.0 DHW thermal excursion calculated using a climatology centered in 1988 is thermally equivalent to a 4.8 DHW excursion calculated using a climatology centered in 2008. Similarly, the 1998 climatological baseline falls halfway between the other two baselines, and the 2008‐equivalent DHW threshold falls halfway between the other two 2008‐equivalent DHW thresholds:

The choice of climatological baseline in the Coral Reef Watch DHW thermal metric is not always made clear, but it is of equal importance to the threshold level (e.g., 4°C‐weeks vs. 8°C‐weeks) in future projections. The above equivalence relationships are derived in the mean over all coral reef locations, and do not capture geographic variations. In this sense they are similar to the DHW threshold framing itself, which already imposes this constraint of global homogeneity.

Statistical Downscaling

We perform statistical downscaling on the coarse‐scale (1°) global model SST projections using the fine‐scale (1 km) MUR SST observational data set. The standard state‐of‐the‐art method for statistical downscaling typically used in ecological projection studies is deterministic, and involves the following simple steps (see, e.g., Van Hooidonk et al. (2016)): (a) At each coarse‐scale model cell, and for each month of the year, estimate the climatology and subtract it from the projected time series, yielding monthly anomaly time series; (b) Interpolate the coarse‐scale monthly anomaly time series onto the fine‐scale (1 km) observational grid; (c) At each fine‐scale pixel, for each month, calculate the climatology using MUR SST data; (d) Add the results of steps 2 and three on a month‐by‐month and pixel‐by‐pixel basis, resulting in fine‐scale projections. This procedure utilizes observational data to construct the fine‐scale climatology and thus can potentially correct systematic bias in the climate model. However, it does not use observations in interpolation (Step 3) but instead assumes deterministic spatial dependence structure across the coarse and fine scales, implying that the coarse‐scale anomalies are downscaled to the fine‐scale grid in a homogeneous way through the time series and spatially. This is a fundamental limitation in the standard downscaling method.

Here, we utilize a novel approach to statistical downscaling, which we describe in greater detail in Ekanayaka et al. (2022). Our motivation was to find a downscaling strategy that had more skill than the standard method described above, and that could produce statistically meaningful uncertainty estimates.

Let y (s ) denote the observational SST at MUR pixel s at month t, for i = 1, …, n, assuming that there are a total n fine‐scale pixels in our study region. Let w (s ) denote the climate model output deterministically interpolated to MUR pixel s , i = 1, …, n. We adopt the statistical downscaling method in Ekanayaka et al. (2022). In particular, we assume: where μ 1, (s ) and μ 2, (s ) represent the large‐scale variation and are modeled as deterministic terms for SST and model output, usually called the trend in geostatistics. Then, we model the joint distribution of {(u 1, (s ), u 2, (s )): i = 1, …, n} by using the basis function representation of a bivariate zero‐mean Gaussian process. In our analysis, we pooled the times series of y (s ) − f (s ) and , where f (s ) represents the output from the standard downscaling procedure, and is the average of interpolated model outputs over the observational years. From these pooled time series, we obtain the empirical orthogonal functions (EOFs). Amongst these functions, we implement the method in Shi and Cressie (2007) and choose EOFs with large absolute‐valued coefficients together with f (s ) and as the trend terms μ 1, (s ) and μ 2, (s ), respectively, but use the remaining to model (u 1, (s ), u 2, (s )) with random coefficients as in Krock et al. (2021). There are several advantages of using such a basis‐function representation: (a) The EOFs in the trend terms are designed to describe systematic spatial departure between observational data and climate model output; (b) The other EOFs with random coefficients enable us to model nonstationary spatial dependence within and between {u 1, (s )} and {u 2, (s )}, thus enabling us to downscale the model output inhomogeneously at different areas (such as coastal regions) in a data‐driven way; (c) Using these basis functions effectively reduces dimensionality and makes our method computationally efficient.

Compared with the standard downscaling method, this novel statistical downscaling method uses observational data in the joint model directly instead of using only their climatology. Our method allows us to simultaneously model the observational data and climate model output, learn their relationship and then use this relationship to produce downscaled projections. Ekanayaka et al. (2022) performed validation studies to compare this method with the standard downscaling method. MUR data before 2018 and climate model output in the Great Barrier Reef region were used as training data to fit the bivariate statistical model. In this methods study performed by our group, we compared the downscaled results from both the standard downscaling method and our new method with withheld “test” MUR data from 2018 to 2020. Over the region containing the entire Great Barrier Reef, we found that the standard downscaling method had mean squared error (MSE) of 0.233°C2 and the BGL method had MSE of 0.214°C2, a reduction of 8%. However, this reduction was more pronounced when averaged only over coral reef locations. Figure 1 presents maps of MSE from the two downscaling methods, in a central region of the Great Barrier Reef. Improvement provided by the BGL downscaling method is especially evident in near‐coastal regions, which is important since many coral reefs globally are located in near‐coastal regions. Averaged over all coral reef locations in this central region including those relatively far from the coast, the standard downscaling method had MSE of 0.252°C2 and the BGL method had MSE of 0.173°C2, a reduction of 31%.

Figure 1

Comparison between standard downscaling and BGL downscaling mean squared error (MSE, in degrees Celsius squared) estimated from validation against withheld 2018–2020 Multiscale Ultrahigh Resolution data in a central region of the Great Barrier Reef. This comparison was performed using SSP126 time series. Coral reef locations are indicated by the brown translucent masking. Note. The MSE improvement provided by the BGL downscaling method that is especially evident in near‐coastal regions. Averaged over coral reef locations, the standard downscaling method had MSE of 0.252°C2 and the BGL method had MSE of 0.173°C2, a reduction of 31%.

BGL also accomplishes our second goal of producing meaningful uncertainty estimates. By using the bivariate statistical model, we are able to quantify the uncertainties associated with the downscaled projections. Note that we obtain from the bivariate model the conditional predictive distribution of y (s )|w (s ) for i = 1, …, n at a future time point t when observational data y (s ) is not available. The downscaled projections are corresponding to the conditional mean, while the conditional standard deviation provides the associated uncertainty. Meanwhile, we note that such uncertainties are based on fitting the model with the training data (i.e., MUR data and climate model output in the observational years) and thus won't be able to characterize uncertainty due to possible extreme departures of the relationship between MUR data and climate model output not presented in the training data in particular unprecedented and unexpected black swan events.

Thermal Departure Projections

We estimate projected times of thermal departure (TD) using the three pairs of DHW thresholds and climatological baselines introduced in Section 2.3. In what follows, we include projections using all three thermal metrics to provide comparability with prior studies, and to quantify the sensitivity of severe bleaching projections to the choice of climatological baseline.

At each 1 km pixel, we concatenate the MUR data from 2002 to 2020 to the mean downscaled projection time series for a particular emissions scenario to create a continuous SST time series from 2002 to 2100. We then calculate the DHW time series from this SST time series, and calculate the year after which every subsequent five year period and every subsequent 10 year period contains at least one heat event surpassing the DHW threshold, at least through 2100. We denote these two TD metrics as TD5Y and TD10Y. Post‐disturbance coral recovery through newly settling recruits requires 7–13 years (Johns et al., 2014) or even >15 years (Baker et al., 2008) if it occurs at all. Thus TD5Y and TD10Y are representative of a range of post‐bleaching coral recovery time scales from damaged but not completely destroyed ecosystems. We note that TD5Y projections might be optimistic, since reefs require more than 5 years to recover after severe bleaching events, but that it is commonly used by prior studies (e.g., Donner (2009); Frieler et al. (2013); Schleussner et al. (2016)). We also note that our construction allows for TD “projections” prior to 2022, and that all TD estimates, even those occurring in the past, depend on information to 2100.

Results

Figure 2 shows the CMIP6 ensemble mean of global mean surface temperature anomaly (GMSTA) over the entire globe in the four emissions scenarios, which begin running in 2014. It also shows the mean of the downscaled SST over all coral reef locations for the four scenarios, including observational MUR data before 2020. Note that the exceptionally strong 2015–2016 El Niño event is clearly apparent in the MUR SST data.

Figure 2

(left) Global mean surface air temperature anomaly projections, relative to an 1880–1900 baseline, from the CMIP6 ensemble mean (right) Mean SST averaged only over coral reef locations included in the analysis, with observational Multiscale Ultrahigh Resolution data before 2020 shown within the shaded region and the downscaled CMIP6 model ensemble projections after 2020. Colors correspond to emissions scenarios as indicated in the legend.

Figure 3 shows global maps for two of the 24 scenarios (4 climate scenarios, 3 DHW metrics, and 2 return timescales) we explored: the highest thermal threshold combination with the latest departure dates and the most optimistic climate scenario (TD5Y, 8 DHW2008, SSP126); and the lowest thermal threshold combination with the earliest departure dates and most pessimistic climate scenario (TD10Y, 8 DHW1988, SSP585). The low‐resolution representations of our high‐resolution results shown in the figures demonstrate general TD dependence on return year, DHW threshold, and cumulative GHG emissions. It is also apparent that some coral reef regions of the world are facing severe thermal stress earlier than others.

Figure 3

Global maps of thermal departure (top) The highest thermal threshold we considered, with the latest departure years, and the most optimistic climate scenario: TD5Y, 8 DHW2008 threshold, and SSP126 (bottom) The lowest thermal threshold we considered, with the earliest departure years, and the most pessimistic climate scenario: TD10Y, 8 DHW1988 threshold, and SSP585. Maps of other scenarios are shown in the Supporting Information S1.

Our main results are shown as cumulative histograms of 1 km2 reef locations remaining under TD5Y and TD10Y (Figure 4) and “slices” through these cumulative histograms at the 30%, 10%, and 1% remaining levels (Tables 1 and 2). Dashes in the tables signify the indicated percent remaining is not crossed before 2100. Vertical gray shading in figures denotes the period of MUR observational data. Note that the drop in reef locations remaining below TD that occurs in ∼2015–2016 corresponds to warming of the reef locations due to the 2015–2016 El Niño visible in the SST data in Figure 2.

Figure 4

Cumulative histograms of thermal departure as a function of year, for SSP126 (black), SSP245 (blue), SSP370 (green), SSP585 (red), for a five year heat event return timescale (TD5Y, top row) and a 10 year heat event return timescale (TD10Y, bottom row). The 1988 (left column) and 2008 (right column) climatological baselines are shown. Cyan and magenta horizontal lines show the 10% and 1% fractional levels respectively; colored vertical ticks on the y‐axis indicate crossings of these levels.

Table 1

Projected Years and GMSTAs After Which Fewer Than the Stated Percentage of 1 km2 Reef Locations Remain Below the Thermal Thresholds, for a Return Timescale of 10 Years (TD10Y)

	8 DHW2008			8 DHW1998			8 DHW1988
	30%	10%	1%	30%	10%	1%	30%	10%	1%
Year in twenty‐first century
SSP126	25	39	__	17	29	__	16	20	34
SSP245	25	35	53	17	28	44	16	18	34
SSP370	26	33	47	19	27	39	16	19	33
SSP585	22	30	42	16	25	36	16	17	30
Global mean surface temperature anomaly (°C)
SSP245	1.4	1.7	1.9	1.2	1.5	1.8	1.1	1.2	1.7
SSP370	1.4	1.7	1.9	1.2	1.5	1.8	1.1	1.2	1.6
SSP585	1.3	1.5	1.9	1.1	1.4	1.7	1.1	1.2	1.5

Table 2

Projected Years and GMSTAs After Which Fewer Than the Stated Percentage of 1 km2 Reef Locations Remain Below the Thermal Thresholds, for a Return Timescale of 5 Years (TD5Y)

	8 DHW2008			8 DHW1998			8 DHW1988
	30%	10%	1%	30%	10%	1%	30%	10%	1%
Year in twenty‐first century
SSP126	30	75	__	23	32	__	19	25	95
SSP245	29	40	62	22	31	49	19	23	40
SSP370	29	36	53	23	30	45	19	25	36
SSP585	26	34	45	21	28	40	19	23	34
Global mean surface temperature anomaly (°C)
SSP245	1.5	1.8	2.0	1.3	1.6	1.9	1.2	1.4	1.8
SSP370	1.5	1.7	2.0	1.4	1.6	1.9	1.2	1.4	1.7
SSP585	1.4	1.6	2.0	1.3	1.5	1.8	1.2	1.4	1.6

It is also useful to interpolate the departure year data using the GMSTA estimates displayed in Figure 2; we perform the interpolation after applying a 10‐year running mean to the GMSTA data. Plots of departure as a function of GMSTA are shown in the Supporting Information S1. Tables 1 and 2 provide GMSTA points of departure beyond various fractions of reefs lost for the four emissions scenarios. Tables 3 and 4 provide percentages and number of reefs remaining below the specified thermal metric, for future GMSTA values.

Table 3

Percentages and Numbers of Reef Locations Remaining Below the Stated Thresholds, for a Return Timescale of 10 Years (TD10Y)

	8 DHW2008			8 DHW1998			8 DHW1988
	1.5°C	1.7°C	2.0°C	1.5°C	1.7°C	2.0°C	1.5°C	1.7°C	2.0°C
Percent 1 km2 reef locations remaining below threshold
SSP245	26%	9%	0%	11%	3%	0%	3%	1%	0%
SSP370	24%	6%	0%	9%	1%	0%	2%	1%	0%
SSP585	15%	3%	0%	5%	1%	0%	1%	0%	0%
Number of 1 km2 reef locations remaining below threshold, out of 773K
SSP245	201K	68 K	4K	83 K	21K	2K	24K	6K	729
SSP370	191K	52 K	9K	73 K	14K	4K	17K	5K	1233
SSP585	117K	25 K	6K	40 K	9K	3K	10K	4K	2265

Table 4

Percentages and Numbers of Reef Locations Remaining Below the Stated Thresholds, for a Return Timescale of 5 Years (TD5Y)

	8 DHW2008			8 DHW1998			8 DHW1988
	1.5°C	1.7°C	2.0°C	1.5°C	1.7°C	2.0°C	1.5°C	1.7°C	2.0°C
Percent 1 km2 reef locations remaining below threshold
SSP245	33%	15%	1%	17%	5%	0%	5%	2%	0%
SSP370	32%	14%	1%	15%	4%	0%	4%	1%	0%
SSP585	21%	6%	1%	9%	2%	0%	2%	1%	0%
Number of 1 km2 reef locations remaining below threshold, out of 773K
SSP245	253K	113 K	7K	132 K	42K	3K	42K	12 K	1250
SSP370	253K	119 K	16K	120 K	36K	6K	34K	11 K	2674
SSP585	171K	50 K	12K	75 K	16K	5K	21K	6K	2628

99% of reef locations are projected to exceed a thermal threshold of 8.0 DHW1988 at least once every 10 years (TD10Y) by 2034, 2034, 2033, and 2030 under SSP126, SSP245, SSP370, and SSP585 (Table 1). In terms of GMSTA, once global heating surpasses 1.5°C–1.7°C, we project that fewer than 1% of reefs will remain below TD10Y, depending on emissions scenario. As of 2021, fewer than 9% of 1 km2 reef locations remained below TD10Y under all emissions scenarios.

TD5Y projections are slightly further in the future than TD10Y projections, as the severe bleaching must occur at least once every five years instead of once every 10 years. 99% of reef locations are projected to exceed TD5Y by 2040, 2036, and 2034 under SSP245, SSP370, and SSP585, corresponding to GMSTAs of 1.8°C, 1.7°C, and 1.6°C, respectively. Higher emissions scenarios push coral reefs over this point at lower GMSTAs due to the progressively steeper rates of global heating (Figure 2), possibly corresponding to less time for deep ocean heat uptake.

As of 2021, fewer than 21% of 1 km2 reef locations remained below TD5Y under all scenarios. We project that at 1.5°C GMSTA, between 2% and 5% of reef locations will remain below TD5Y, and between 1% and 3% will remain below TD10Y. We project that at 2.0°C GMSTA, the number of reef locations remaining below TD5Y or TD10Y (fewer than 2700 and 2300 1 km2 locations respectively) will be closer to 0% than to 1%.

Under all the thermal metrics, the SSP126 scenario, although still dire, projects a markedly better prognosis for corals than the other three emissions scenarios. Under TD5Y, 1% of reefs are projected to remain below the thermal threshold until 2095. Also, although 99% of reefs surpass the threshold under TD10Y by 2034, further losses proceed more slowly than in the other three emissions scenarios (Figure 4).

We validated our analysis by comparing the mean of the three annual maximum ocean heat events at each reef pixel from 2018 to 2020 in the downscaled SSP126 SST time series to the corresponding value in the MUR SST data. We found that the mean of a distribution of MUR values subtracted from corresponding downscaled model SST values was −1.8°C‐weeks (with a standard deviation of 1.7°C‐weeks), that is, the downscaled model value underestimated the MUR data by 1.8°C‐weeks (see Figure S7 in Supporting Information S1). We found similar results for the other three SSPs. This suggests that the projections are “conservative” in the sense that they underestimate future coral bleaching.

Discussion and Conclusion

In 2020, global heating (GMSTA) was 1.2°C–1.3°C above pre‐industrial levels, and human GHG emissions will likely push Earth to 1.5°C GMSTA sometime in the 2030s, according to CMIP6 model projections (Figure 2). Unless humanity accomplishes climate mitigation approximating the SSP126 scenario, Earth will likely surpass 2°C GMSTA around mid‐century (e.g., Table 1). We have provided projections, with unprecedented spatial resolution, of future years and global heating levels beyond which coral severe bleaching conditions due to this anthropogenic global heating will be continuous relative to coral recovery timescales. Novel aspects of our departure year and GMSTA projections include using the CMIP6 model ensemble; attaining 1 km resolution; downscaling with an improved method; performing an end‐to‐end validation against observational data; and providing projections under six combinations of two ecologically relevant severe bleaching event return timescales (5 and 10 years) and three DHW thresholds.

Clarifying that complete specification of DHW thresholds requires not one, but two numbers facilitates apples‐to‐apples comparisons with prior studies. Schleussner et al. (2016) projected a 70%–90% loss at 1.5°C and 99% loss at 2°C GMSTA, using CMIP3 global models (without downscaling) and a thermal criteria of TD5Y and 8 DHW1990 (the center of a 1980–2000 reference climatology). These results were adopted by the IPCC Special Report on Global Warming of 1.5°C (Summary for policymakers, 2018). Using nearly identical thermal criteria (TD5Y and 8 DHW1988), we project a 95%–98% loss at 1.5°C and a 99.7% loss at 2°C GMSTA (Table 4).

Donner (2009) used one global model and a thermal metric of TD5Y and 8 DHW1988 (a 1985–2000 climatology) to project roughly 70% of coarse‐scale (not downscaled) global model locations will surpass the metric in 2025%, and 90% by 2040, under SRES B1 (similar to SSP245); our study projects 2019 and 2023 (Table 2).

Frieler et al. (2013), using 19 CMIP3 models and an 8 DHW1990 (1980–1999 climatology), found that 90% of coarse grid cells surpass TD5Y at 1.5°C, and that all grid cells surpass TD5Y before 2°C GMSTA; our study projects over 95% TD5Y at 8 DHW1988 and 1.5°C, and over 99.7% at 2°C (Table 4).

Van Hooidonk et al. (2016) was the only prior study that applied statistical downscaling; they downscaled CMIP5 projections to 4 km resolution and found mean TD1Y values (annual recurrence) of ocean heat events surpassing 8 DHW1995 (1982–2008 climatology) of 2054 for the climate scenarios RCP 4.5 and 2043 for RCP 8.5, which are similar to the scenarios SSP245 and SSP585 used here. Our study does not include comparable metrics, and we note that annual severe bleaching might be too “conservative” a metric to be useful, given observed post‐bleaching recovery times of about a decade.

Dixon et al. (2022) applied statistical downscaling with MUR data to 1 km resolution and used the low thermal criteria of TD10Y and 4 DHW1988 and found 0.2% of reefs at 1.5°C, and 0% at 2.0°C of mean global heating. We do not use such a low thermal threshold, so we cannot perform a comparison, but we note again that projections of remaining reefs depend critically on choice of thermal threshold.

Our results project an earlier decline for the world's coral reefs than either Schleussner et al. (2016) or Donner (2009), but are in agreement with Frieler et al. (2013). However, these earlier studies used a 5‐year return timescale, but a 10‐year return timescale is more ecologically appropriate.

There are three realms of uncertainty in our projections. The first is scenario uncertainty, the uncertainty over humanity's collective future emissions; this dimension is spanned over the four “SSP” emissions scenarios. The second realm of uncertainty is projection uncertainty, part of which stems from uncertainties in the global models (Lehner et al., 2020). Projection uncertainty, in the context of ecological projections, can also arise from uncertainties in observational datasets and from the downscaling methodology. The two prior studies that do estimate projection uncertainty do so from the spread of individual global models within the model ensemble (Frieler et al., 2013; Schleussner et al., 2016). However, we cannot apply this method directly to our downscaled results. One key area for future work is to understand and reduce projection uncertainty. We are currently developing a statistical uncertainty quantification from the BGL downscaling method and the model ensemble (informed by comparative assessments between individual models and observations). In addition to uncertainty quantification, skill‐weighting the ensemble could allow better use of information, potentially improving projection accuracy, which could be checked in hindcast experiments. Furthermore, the current standard practice of using what amounts to an arbitrary collection of models and taking their ensemble means creates uncertainty. To illustrate this, we performed our analysis on a separate CMIP6 ensemble of 127 model members from 27 model groups (Supporting Information Text T2 and Tables S1 and S2 in Supporting Information S1). The different ensemble led to slightly different results, for example, projecting 2% of reef locations to not surpass 8 DHW1988 at TD10Y under SSP245, as opposed to 3%. The 127‐member ensemble projects 99% of reefs to exceed 8 DHW1988 at TD10Y under SSP126 in 2086, as compared to 2034 for the 35‐member ensemble; this seemingly dramatic difference can be explained by the flattening of the cumulative histogram curve in bottom left panel of Figure 4 near this 99% remaining threshold, creating a large rate of change. More serious is the possibility of misidentifying specific locations of projected refugia. This arbitrariness due to the ensemble composition, which is also present in other ensemble studies, could be eliminated in the future via skill‐weighting.

The third realm of uncertainty is ecological uncertainty, the uncertainty in the relationship between ocean heat events and the response of coral reefs. We have spanned a small part of this realm by providing projections under the two severe bleaching recovery timescales, and three thermal threshold metrics.

As is the case with the prior studies, our study does not factor in additional ecologic factors which could potentially mitigate or exacerbate coral reef degradation and loss. On shorter timescales, clouds can block sunlight, potentially reducing algal production of reactive oxygen species (M. E. Baird et al., 2018; Skirving et al., 2018; Roth, 2014), and mitigating bleaching during marine heat events (Mumby et al., 2001). Reef depth could also affect bleaching by reducing sunlight and water temperatures (A. H. Baird et al., 2018; Frade et al., 2018; Muir et al., 2017; Smith et al., 2014). Relatively high SST variability correlates with lower bleaching risk (Beyer et al., 2018; Safaie et al., 2018). Relatively high nutrient levels correlates with higher bleaching risk (DeCarlo & Harrison, 2019).

On longer timescales, dispersal of coral larvae could result in establishment of populations in cooler regions of the future ocean (Greenstein & Pandolfi, 2008). Ocean acidification, sea‐level‐rise, sedimentation, and intensifying storms could further harm corals (Blanchon et al., 2009; Cheal et al., 2017; Cohen et al., 2009; Field et al., 2011; Hoegh‐Guldberg et al., 2007; Perry et al., 2018).

In this study, we do not attempt to account explicitly for highly uncertain coral adaptation, although our use of three climatological baselines could serve as a rudimentary proxy. Adaptation of corals and/or symbionts (such as acclimatization, symbiont shuffling, or genetic change) would improve coral prospects, but evidence is equivocal and mechanisms remain poorly understood (Baker et al., 2004; Chakravarti et al., 2017; Donner et al., 2005; Hoegh‐Guldberg, 2014; Parmesan, 2006; Torda et al., 2017). Logan et al. (2021) folds potential symbiont‐mediated adaptive capacity from symbiont shuffling and symbiont evolution into thermal viability projections from an ecological model, driven by SST output from a global climate model. Shuffling of symbionts with assumed thermal growth optima of up to 1.5°C above heat‐sensitive symbionts allowed the model to simulate thriving global reefs beyond 2100. Even under the most extreme climate scenario (RCP 8.5), 23% of simulated global reefs remained healthy under symbiont shuffling combined with symbiont evolution.

A major focus for future work will be understanding and constraining ecological uncertainty. Adaptation can be included in coral projections when based on observed adaptation levels, as hypothetical adaptation levels lead to unconstrained projections. It might also be possible to constrain the coral response to ocean heat events through the use of empirical data, such as remotely sensed severe coral bleaching from satellite platforms. This could provide sufficient data to create models of the coral response that account for the coral locations, and could include additional predictor variables.

Our analysis does provide projected 1 km2 locations of global coral refugia. However, given the high degree of uncertainty, and imminent data science innovations with the potential to constrain this uncertainty, we choose not to highlight the identification of refugia in our current study, despite having created an online visualizer. We note that a small number of reefs are projected to persist beyond 2°C GMSTA even under the most stringent metric (Table 3), but that we have low confidence in the precise locations of these potential refugia. Indeed, we see an urgent need to further improve ecological projection in order to attain the capacity to robustly identify refugia, including understanding the physical basis for their projected persistence, for the sake of guiding conservation efforts. Our group plans to release improved projections in a subsequent study, which will include identification of refugia.

Finally, we feel that it is no longer possible to overstate the importance of rapid cessation of human GHG emissions. In the absence of extremely rapid coral adaptation to increasing heat, which would need to occur in the simultaneous presence of the many additional and serious anthropogenic stressors listed earlier, our results suggest that 2°C of global heating could render Earth essentially uninhabitable to warm water coral reefs as we know them. Furthermore, if near‐future emissions are equivalent or greater than SSP245, we project that by 2040 over 99% of the world's reefs will be subject to thermal severe bleaching conditions too recurrent for recovery (TD5Y), which will continue to worsen. On the other hand, if emissions approximated the SSP126 scenario and GMSTA were limited to 1.5°C, this level of severe bleaching might not attain and global conditions could stabilize on a planet with coral reefs.

Supporting Information S1

Click here for additional data file.