# Libraries
  library(lubridate)
  library(tidyverse)
  library(gridExtra) 
  library(zoo)
  library(reshape2)
  library(knitr)

# Default ggplot settings
  ggthe_bw<-theme(plot.background=element_blank(),
          panel.grid.major.y = element_blank(),
          panel.grid.major.x = element_blank(),
          panel.grid.minor.x = element_blank(),
          panel.grid.minor.y = element_blank(),
          legend.box.background = element_rect(),
          panel.background =element_rect(fill = NA, color = "black")
          )+
    theme_bw()

1. Temperature data sources

Obtain temperature data from the different sources

#OI_SST (provided by Ruben van Hooidonk)
    Uva.OI<-read.csv("B.Temperature_data/Uva_Daily_OISST_1982-2016.csv",header = T)
    Uva.OI$Date<-as.Date(Uva.OI$Date, format="%Y-%m-%d")

#InSitu (provided by Tyler Smith)
    Uva.InSitu<-read.csv("B.Temperature_data/Uva_Daily_InSitu_1997-2018.csv", header = T)
    Uva.InSitu$Date<-as.Date(Uva.InSitu$Date, format="%Y-%m-%d")

# CRW-5km SST (provided by Ruben van Hooidonk)
    Uva.CRW<-read.csv("B.Temperature_data/Uva_Daily_CRW_5km_1985-2016.csv",header = T)
    Uva.CRW$Date<-as.Date(Uva.CRW$Date, format="%Y-%m-%d")
    
#CCI_SST (provided by Ruben van Hooidonk)
    Uva.CCI<-read.csv("B.Temperature_data/Uva_Daily_CCI_1982-2016.csv",header = T)
    Uva.CCI$Date<-as.Date(Uva.CCI$Date, format="%Y-%m-%d")

Compile all temperature data in wide format

Temperature.Wide<-merge(Uva.OI, Uva.InSitu, on=Date, all=TRUE)
Temperature.Wide<-merge(Temperature.Wide, Uva.CRW, on=Date, all=TRUE)
Temperature.Wide<-merge(Temperature.Wide, Uva.CCI, on=Date, all=TRUE)

Compile all temperature data in tall format

#CCI_SST
    Uva.CCI.L<-Uva.CCI
    Uva.CCI.L$Source<-"CCI_SST"
    colnames(Uva.CCI.L) <- c("Date","Temperature","Source")

#OISST
    Uva.OI.L<-Uva.OI
    Uva.OI.L$Source<-"OI_SST"
    colnames(Uva.OI.L) <- c("Date","Temperature","Source")

#InSitu
    Uva.InSitu.L<-Uva.InSitu
    Uva.InSitu.L$source<-"In_situ"
    colnames(Uva.InSitu.L) <- c("Date","Temperature", "Source")
    
# CRW-5km SST
      Uva.CRW.L<-Uva.CRW
      Uva.CRW.L$Source<-"CRW_SST"
      colnames(Uva.CRW.L) <- c("Date", "Temperature", "Source")
      
Temperature.Tall<-rbind(Uva.OI.L, Uva.InSitu.L, Uva.CRW.L, Uva.CCI.L)

2. Climatology and MMM

2.1 OI MMM

Calculate the mean SST per month-year, the monthly mean and then obtain the maximum monthly mean (MMM) using data OISST from 1985-2012.

# Extract months and years from date
  Uva.OI$Month<-month(Uva.OI$Date)
  Uva.OI$Year<-year(Uva.OI$Date)

# Select data from 1985-2012
  Uva.OI.1985_2012<-filter(Uva.OI, Year>1984)
  Uva.OI.1985_2012<-filter(Uva.OI.1985_2012, Year<2013)

# Calculate month-year mean SST
  Uva.OI.MM_1985_2012<-Uva.OI.1985_2012 %>%
    group_by(Month, Year) %>%
    summarize(MM = mean(OISST, na.rm = TRUE))

# Calculate monthly mean SST
  Uva.OI.MonthlyClimatology_1985_2012<-Uva.OI.MM_1985_2012 %>%
    group_by(Month) %>%
    summarize(Clima_M = mean(MM, na.rm = TRUE))
  Uva.OI.MonthlyClimatology_1985_2012
# Get maximum monthly mean (MMM)
  MMM_OI<-max(Uva.OI.MonthlyClimatology_1985_2012$Clima_M)
  MMM_OI
## [1] 29.16445

2.2 CRW MMM

Import the MMM from CRW 5km suit products

# Calculated MMM
# # Extract months and years from date
#   Uva.CRW$Year<-year(Uva.CRW$Date)   
#   Uva.CRW$Month<-month(Uva.CRW$Date)   
# 
# # Calculate month-year mean SST
#   Uva.CRW.MM<-Uva.CRW %>%
#     group_by(Month, Year) %>%
#     summarize(M_Temeprature = mean(Temperature, na.rm = TRUE))
# 
# # Calculate monthly mean SST             
#   Uva.CRW.MonthlyClimatology<-Uva.CRW.MM %>%
#     group_by(Month) %>%
#     summarize(Clima_M = mean(M_Temeprature, na.rm = TRUE))
#   Uva.CRW.MonthlyClimatology
#   
# # Get maximum monthly mean (MMM)              
#   MMM_CRW<-max(Uva.CRW.MonthlyClimatology$Clima_M)
#   MMM_CRW

# Get maximum monthly mean (MMM) from CRW 5km products
  Uva.CRW.MMM<-read.csv("B.Temperature_data/Uva_CRW_MMM_1985-2012.csv", header = T)
  
  MMM_CRW<-Uva.CRW.MMM[1,1]
  MMM_CRW
## [1] 28.9

2.3 In Situ MMM

This code uses in situ available data collected by Tyler Smith et al and calculates an approximated In Situ climatology (~1997-2018, with gaps). The time frame used in this climatology differs from the CWR and OI climatologies. We decided to use all the data available since data were scarcer.

# Extract months and years from date
  Uva.InSitu$Year<-year(Uva.InSitu$Date)   
  Uva.InSitu$Month<-month(Uva.InSitu$Date)   

# Calculate month-year mean SST
  Uva.InSitu.MM<-Uva.InSitu %>%
    group_by(Month, Year) %>%
    summarize(M_Temeprature = mean(Temperature, na.rm = TRUE))

# Calculate monthly mean SST             
  Uva.InSitu.MonthlyClimatology<-Uva.InSitu.MM %>%
    group_by(Month) %>%
    summarize(Clima_M = mean(M_Temeprature, na.rm = TRUE))
  Uva.InSitu.MonthlyClimatology
# Get maximum monthly mean (MMM)              
  MMM_InSitu<-max(Uva.InSitu.MonthlyClimatology$Clima_M)
  MMM_InSitu
## [1] 29.05983

2.4 CCI MMM

Calculate the mean SST per month-year, the monthly mean and then obtain the maximum monthly mean (MMM) using data CCI from 1985-2012.

# Extract months and years from date
  Uva.CCI$Month<-month(Uva.CCI$Date)
  Uva.CCI$Year<-year(Uva.CCI$Date)

# Select data from 1985-2012
  Uva.CCI.1985_2012<-filter(Uva.CCI, Year>1984)
  Uva.CCI.1985_2012<-filter(Uva.CCI.1985_2012, Year<2013)

# Calculate month-year mean SST
  Uva.CCI.MM_1985_2012<-Uva.CCI.1985_2012 %>%
    group_by(Month, Year) %>%
    summarize(MM = mean(SST_CCI, na.rm = TRUE))

# Calculate monthly mean SST
  Uva.CCI.MonthlyClimatology_1985_2012<-Uva.CCI.MM_1985_2012 %>%
    group_by(Month) %>%
    summarize(Clima_M = mean(MM, na.rm = TRUE))
  Uva.CCI.MonthlyClimatology_1985_2012
# Get maximum monthly mean (MMM)
  MMM_CCI<-max(Uva.CCI.MonthlyClimatology_1985_2012$Clima_M)
  MMM_CCI
## [1] 29.24618

Notes about MMMs

  • Values: CRW MMM (28.9 C) < Insitu MMM (29.05983 C) < OI MMM (29.16445 C) < CCI (29.24618)
  • Observations: Insitu MMM (n=7,864 days) < CRW MMM (n=10,227) = OI MMM (n=10,227) = CCI MMM (n=10,227)
  • Time frame: Insitu MMM (1997-2018, with gaps) < CRW MMM (1985-2012) = OI MMM (1985-2012) = OI MMM (1985-2012)

Calculated CRW MMM is higher (29.09882 C) than MMM extracted directly from the 5km products (28.9 C), (ftp://ftp.star.nesdis.noaa.gov/pub/sod/mecb/crw/data/5km/v3.1/climatology/nc/ct5km_climatology_v3.1.nc).

https://coralreefwatch.noaa.gov/satellite/bleaching5km/index.php says: This 28-year climatology covers the time period 1985-2012 and is based on the daily global 5km ‘CoralTemp’ SST data product. This 5km climatology was re-centered to the time-center of the operational 50km climatology (1985-1990 plus 1993) to maintain consistency with the former Version 2 5km climatology and CRW’s heritage 50km products. Not really sure how the data are “re-centered”, but this seems to lower the MMM.

3. DHW

3.1 OISST DHW calculation

This code uses a mixed approach (Ruben van Hooidonk + CRW) to estimate DHW, based on OI_SST data:

  • SST data source is OI_SST v2, not CRW (CoralTemp). OI data were initially preferred because CRW data start in 1985 and missed the 1982-83 ENSO
  • MMM values used are the calculated above, based on OI_SST_v2 SST from the period 1985-2012. This is the same period used for the CRW_5KM_V3_climatology.
  • OISST has one temperature value per day, there is no distinction between day/night SST as there is in CRW products
  • HotSpots are calculated daily, then divided by 7 to calculate DHW from DHDays
##### Calculate IO DHW
      Uva.OI.L$HotSpot<-(Uva.OI.L$Temperature-(MMM_OI))
      Uva.OI.L$HotSpot<-ifelse(Uva.OI.L$HotSpot>=0,
                              Uva.OI.L$HotSpot, 0)# REMOVE NEGATIVE ANOMALIES
      Uva.OI.L$Stress<-ifelse(Uva.OI.L$HotSpot>=1,
                              Uva.OI.L$HotSpot, 0) # REMOVE HotSpot lower than 1
      Uva.OI.L$W_Stress<-(Uva.OI.L$Stress/7) # DHDays to DHWeeks
      Uva.OI.L$DHW<-rollapplyr(Uva.OI.L$W_Stress,list(-(83:0)),sum,fill=NA)# Add 12 weeks of heat stress data

3.2 CRW SST DHW calculation

This code uses CRW MMM and CRW SST data to estimate DHW:

  • MMM values used were extracted from the CWR 5km climatology that includes 1985-2012 data. Calculated MMM value is slightly different (warmer) than MMM value extracted from CRW_5KM_V3_climatology.
##### Calculate CRW DHW based on CRW MMM
      Uva.CRW.L$HotSpot<-(Uva.CRW.L$Temperature-(MMM_CRW))
      Uva.CRW.L$HotSpot<-ifelse(Uva.CRW.L$HotSpot>=0,
                              Uva.CRW.L$HotSpot, 0)# REMOVE NEGATIVE ANOMALIES
      Uva.CRW.L$Stress<-ifelse(Uva.CRW.L$HotSpot>=1,
                              Uva.CRW.L$HotSpot, 0) # REMOVE HotSpot lower than 1
      Uva.CRW.L$W_Stress<-(Uva.CRW.L$Stress/7) # DHDays to DHWeeks
      Uva.CRW.L$DHW<-rollapplyr(Uva.CRW.L$W_Stress,list(-(83:0)),sum,fill=NA)# Add 12 weeks of heat stress data

3.3 InSitu DHW calculation

This code uses a mixed approach (Ruben van Hooidonk + CRW) to estimate DHW, based on In situ data:

  • Temperature data source are Tyler Smith’s Hobo temps, not CRW SST.
  • MMM values used are the calculated above, based HOBOs’ data from the period ~ 1997-2018. This is a different period than the one used for the CRW_5KM_V3_climatology, but it is the range available. The data present gaps
  • Daily mean temperatures were calculated using all the data available for a day (day + night)
  • HotSpots are calculated daily, then divided by 7 to calculate DHW from DHDays
##### Calculate IO DHW
      Uva.InSitu.L$HotSpot<-(Uva.InSitu.L$Temperature-(MMM_InSitu))
      Uva.InSitu.L$HotSpot<-ifelse(Uva.InSitu.L$HotSpot>=0,
                              Uva.InSitu.L$HotSpot, 0)# REMOVE NEGATIVE ANOMALIES
      Uva.InSitu.L$Stress<-ifelse(Uva.InSitu.L$HotSpot>=1, 
                              Uva.InSitu.L$HotSpot, 0) # Select HotSpots >=1
      Uva.InSitu.L$W_DHW<-(Uva.InSitu.L$Stress/7)
      Uva.InSitu.L$DHW<-rollapplyr(Uva.InSitu.L$W_DHW,list(-(83:0)),sum,fill=NA)

3.4 CCI DHW calculation

This code uses a mixed approach (Ruben van Hooidonk + CRW) to estimate DHW, based on CCI_SST data:

  • SST data source is ESA Sea Surface Temperature Climate Change Initiative (ESA SST CCI): GHRSST Multi-Product ensemble (GMPE) ftp://anon-ftp.ceda.ac.uk/neodc/esacci/sst/data/CDR_v2/Analysis/L4/v2.1.
  • CCI data were added since someone told us the OISST is not so good on early years. So we need to test the effect on the DHW calculations.
  • MMM values used are the calculated above, based on CCI_v2.1 SST from the period 1985-2012. This is the same period used for the CRW_5KM_V3_climatology.
  • Ruben gave me one temperature value per day (mean?), there is no distinction between day/night SST as there is in CRW products
  • HotSpots are calculated daily, then divided by 7 to calculate DHW from DHDays
##### Calculate IO DHW
      Uva.CCI.L$HotSpot<-(Uva.CCI.L$Temperature-(MMM_CCI))
      Uva.CCI.L$HotSpot<-ifelse(Uva.CCI.L$HotSpot>=0,
                              Uva.CCI.L$HotSpot, 0)# REMOVE NEGATIVE ANOMALIES
      Uva.CCI.L$Stress<-ifelse(Uva.CCI.L$HotSpot>=1,
                              Uva.CCI.L$HotSpot, 0) # REMOVE HotSpot lower than 1
      Uva.CCI.L$W_Stress<-(Uva.CCI.L$Stress/7) # DHDays to DHWeeks
      Uva.CCI.L$DHW<-rollapplyr(Uva.CCI.L$W_Stress,list(-(83:0)),sum,fill=NA)# Add 12 weeks of heat stress data

4. DHW plots

4.1 Explore available data based on all data sources

# Wide DHW 
    DHW.InSitu.W<-Uva.InSitu.L[, c("Date","DHW")]
    colnames(DHW.InSitu.W) <- c("Date","InSitu_DHW")
    
    DHW.CRW.W<-Uva.CRW.L[, c("Date","DHW")]
    colnames(DHW.CRW.W) <- c("Date","CRW_DHW")
    
    DHW.OI.W<-Uva.OI.L[, c("Date", "DHW")]
    colnames(DHW.OI.W) <- c("Date","OI_DHW")
    
    DHW.CCI.W<-Uva.CCI.L[, c("Date", "DHW")]
    colnames(DHW.CCI.W) <- c("Date","CCI_DHW")
    
    DHW.Wide<-merge(DHW.InSitu.W, DHW.CRW.W, on=Date, all=TRUE)
    DHW.Wide<-merge(DHW.Wide, DHW.OI.W, on=Date, all=TRUE)
    DHW.Wide<-merge(DHW.Wide, DHW.CCI.W, on=Date, all=TRUE)
  
# Long DHW data
    Uva.InSitu.DHW<-Uva.InSitu.L[, c("Source", "Date", "Temperature", "HotSpot", "DHW")]
    Uva.CRW.DHW<-Uva.CRW.L[, c("Source", "Date", "Temperature",  "HotSpot", "DHW")]
    Uva.OI.DHW<-Uva.OI.L[, c("Source", "Date", "Temperature",  "HotSpot", "DHW")]
    Uva.CCI.DHW<-Uva.CCI.L[, c("Source", "Date", "Temperature",  "HotSpot", "DHW")]

    Long.data.DHW<-rbind(Uva.InSitu.DHW, Uva.CRW.DHW, Uva.OI.DHW, Uva.CCI.DHW)
    Long.data.DHW$Source<-factor(Long.data.DHW$Source, 
                            levels = c("CCI_SST", "OI_SST","CRW_SST","In_situ"))
    Long.data.DHW$Year<-lubridate::year(Long.data.DHW$Date)
    Long.data.DHW$Month<-lubridate::month(Long.data.DHW$Date)

# Highest DWH per year
    MaxDHW<-as.data.frame(Long.data.DHW %>%
          group_by(Source, Year) %>%
          summarize(max = max(DHW, na.rm = TRUE)))
    MaxDHW<-dcast(MaxDHW, Year ~ Source, value.var = "max")
    MaxDHW
# Highest HotSpot per year
  MaxHS<-as.data.frame(Long.data.DHW %>%
        group_by(Source, Year) %>%
        summarize(max = max(HotSpot, na.rm = TRUE)))
  MaxHS<-dcast(MaxHS, Year ~ Source, value.var = "max")
  MaxHS
# Plot data
AllSources<-ggplot(Long.data.DHW, aes(x=Date)) + 
    ggthe_bw+ theme(legend.position = "none")+
    facet_wrap(~Source)+
    geom_hline(yintercept=8, linetype="dashed", color="red")+
    geom_hline(yintercept=4, linetype="dashed", color="orange")+
    geom_line(aes(y = Temperature, colour = "Temperature")) +
    geom_line(aes(y = DHW, colour = "DWH")) +
    ylab("DWH (C-week)  /  Temperature (C)")
AllSources

4.2 Comparison of the 3 ENSO events based on OISST data

# 1. Extract only 1982-83, 1997-98 and 2015-16 Data 

  # 1981-82
    Days<-as.data.frame(seq(1,730,by=1))
    names(Days) <- c("Days")
    
    N_1982<-as.data.frame(seq(as.Date("1982-01-01"), 
                              as.Date("1983-12-31"), by="days"))
        names(N_1982) <- c("Date")
        N_1982<-cbind(Days, N_1982)
        N_1982$Event<-"1982-83"
        tail(N_1982,n=2)
  # 1997-98
    N_1997<-as.data.frame(seq(as.Date("1997-01-01"), 
                              as.Date("1998-12-31"), by="days"))
        names(N_1997) <- c("Date")
        N_1997<-cbind(Days,N_1997)
        N_1997$Event<-"1997-98"
        tail(N_1997,n=2)
  # 2015-16  
    Days<-as.data.frame(seq(1,731,by=1)) # 2015-2016 Needs an extra day 
    names(Days) <- c("Days")
    
    N_2015<-as.data.frame(seq(as.Date("2015-01-01"), 
                              as.Date("2016-12-31"), by="days"))
        names(N_2015) <- c("Date")
        N_2015<-cbind(Days, N_2015)
        N_2015$Event<-"2015-16"
        tail(N_2015,n=2)
  # ENSO years
    ENSO_Dates<-rbind(N_1982, N_1997, N_2015)

Nino.OI_data<-merge(x = ENSO_Dates, y = Uva.OI.L, by = "Date", all.x = TRUE)
Nino.OI_data$Dates<-format(Nino.OI_data$Date, format= "%m-%d")

Nino.CCI_data<-merge(x = ENSO_Dates, y = Uva.CCI.L, by = "Date", all.x = TRUE)
Nino.CCI_data$Dates<-format(Nino.CCI_data$Date, format= "%m-%d")
#head(Nino.OI_data,n=2)
#tail(Nino.OI_data,n=2)

Figure 1a: 3 ENSOs

Together

Figure1a<-ggplot(Nino.OI_data, aes(x=Days)) + ggtitle("a.")+ ggthe_bw+
  theme(legend.box.background = element_rect(),
        legend.position = c(0.15, 0.75)) +
  geom_hline(yintercept=8, linetype="dashed", color="red", size=1)+
  geom_hline(yintercept=4, linetype="dashed", color="orange", size=1)+


  geom_line(aes(y = DHW), size=1.0) + 
  scale_x_continuous(limits = c(0,701), expand = c(0, 0), 
                     name=("Month"),
                     breaks = seq(0,701, by= 90))+
  scale_y_continuous(limits = c(0,11), 
                     name=(expression("DHW"~(degree*C ~"-weeks"))), expand = c(0.01, 0.01),
                     breaks = seq(0,10, by=2)) +  
  facet_grid(~Event)
Figure1a

#ggsave(file="C.Outputs/Figure_1a2.svg", plot=Figure1a2, width=7, height=4)

Figure 1a: Accumulated heat stress and thermal regime at Uva Island, Gulf of Chiriquí. a. Degree heating weeks (DHW) during the 1982-83, 1997-98 and 2015-16 ENSOs estimated using daily OISST data. The yellow dashed line marks the 4 DHW threshold for significant likelihood of bleaching, and the red dashed line marks 8 DHW threshold for significant likelihood of mortality.

4.3 Temperature and DHW during 2015-16 ENSO events based on all data sources OI and InSitu data

Figure 3a

 # Figure 3a

TemperaturePlot<-ggplot(Temperature.Wide, aes(x=Date)) + ggthe_bw+
    theme(legend.position = "none") +
    scale_colour_manual(values = c("gray", "black"), name="Source")
    #scale_colour_manual(values = c("black", "blue4","grey"), name="Source")

TemperaturePlot<-TemperaturePlot + 
  geom_hline(yintercept=MMM_OI+1, linetype="dashed", color="black")+
  geom_hline(yintercept=MMM_InSitu+1, linetype="dashed", color="grey")+
  #geom_hline(yintercept=MMM_CRW+1, linetype="dashed", color="blue4")+
  geom_line(aes(y = (OISST), colour = "OI_SST"), linetype=1, size=0.3, alpha=0.9) +      
  geom_line(aes(y = (Temperature), colour = "In situ"), size=0.4, alpha=0.9) +
  #geom_line(aes(y = (CRWSST), colour = "CRW_SST"), linetype=1, size=0.3, alpha=0.9) +
  
  scale_x_date(date_breaks = "3 months", date_labels = "%y-%m", 
               limits = c(as.Date("2014-07-15"), as.Date("2016-07-10")),
               name = "") +
  scale_y_continuous(limits = c(26,31.1), expand = c(0.01, 0), 
                   name=(expression("Temperature"~(degree*C))),
                   breaks = seq(26,31, by=1))+
  annotate("rect", xmin= as.Date("2014-08-17"), xmax = as.Date("2014-08-24"), 
                 ymin = 26, ymax = 31, fill="blue", alpha = .3) +
  annotate("rect", xmin= as.Date("2015-08-16"), xmax = as.Date("2015-08-22"), 
                 ymin = 26, ymax = 31, fill="blue", alpha = .3) + 
  annotate("rect", xmin= as.Date("2016-04-16"), xmax = as.Date("2016-04-21"), 
                 ymin = 26, ymax = 31, fill="blue", alpha = .3)
  
DHWPlot<-ggplot(DHW.Wide, aes(x=Date)) + ggthe_bw+
   theme(legend.box.background = element_rect(),
        legend.position = c(0.2, 0.4)) +
  scale_colour_manual(values = c("gray", "black"), name="Data source")
  #scale_colour_manual(values = c("blue4", "gray","black"), name="Source")
DHWPlot<-DHWPlot +
    geom_hline(yintercept=8, linetype="dashed", color="red", size=0.5)+
    geom_hline(yintercept=4, linetype="dashed", color="orange", size=0.5)+
    geom_line(aes(y = (OI_DHW), colour = "OI SST"), linetype=1, size=0.3, alpha=0.9) +      
    geom_line(aes(y = (InSitu_DHW), colour = "In situ"), size=0.4, alpha=0.9) +
    #geom_line(aes(y = (CRW_DHW), colour = "CRW"), size=0.4, alpha=0.9) +  
    scale_x_date(date_breaks = "3 months", date_labels = "%m/%Y", 
               limits = c(as.Date("2014-07-15"), as.Date("2016-07-10")),
               name = "Month-Year") +
    scale_y_continuous(limits = c(0,10), expand = c(0.01, 0), 
                   name=(expression("DHW"~(degree*C~"-weeks"))),
                   breaks = seq(0,10, by=2))+
    annotate("rect", xmin= as.Date("2014-08-17"), xmax = as.Date("2014-08-24"), 
                 ymin = 0, ymax = 10, fill="blue", alpha = .3) +
    annotate("rect", xmin= as.Date("2015-08-16"), xmax = as.Date("2015-08-22"), 
                 ymin = 0, ymax = 10, fill="blue", alpha = .3) + 
    annotate("rect", xmin= as.Date("2016-04-16"), xmax = as.Date("2016-04-21"), 
                 ymin = 0, ymax = 10, fill="blue", alpha = .3)

Figure2a_Final<-grid.arrange(TemperaturePlot, DHWPlot, nrow=2, heights=c(6/10, 4/10))

#ggsave(file="C.Outputs/Figure_1b.svg", plot=Figure2a_Final, width=6.5, height=4.5)

Figure 3a: Accumulated heat stress and thermal regime at Uva Island, Gulf of Chiriquí. a. Thermal regime before and during the 2015-16 El Niño using daily OISST (black) and in situ (gray) data collected at Uva Island reef (~3m depth). Upper plot and right-hand axis show mean daily temperature (solid lines) and bleaching threshold (dashed lines) for each data set. Lower plot and left-hand axis show DHW calculated from each data source. The blue vertical lines represent the periods when the coral samples were collected. 

 # Figure supplementary

TemperaturePlotAll<-ggplot(Temperature.Wide, aes(x=Date)) + ggthe_bw+
    theme(legend.position = "none") +
    scale_colour_manual(values = c("gray", "black", "purple", "red"), name="Source")
   
TemperaturePlotAll<-TemperaturePlotAll + 
  geom_hline(yintercept=MMM_OI+1, linetype="dashed", color="purple")+
  geom_hline(yintercept=MMM_InSitu+1, linetype="dashed", color="black")+
  geom_hline(yintercept=MMM_CRW+1, linetype="dashed", color="grey")+
  geom_hline(yintercept=MMM_CCI+1, linetype="dashed", color="red")+
  
  geom_line(aes(y = (OISST), colour = "OI_SST"), linetype=1, size=0.3, alpha=0.9) +      
  geom_line(aes(y = (Temperature), colour = "In situ"), size=0.4, alpha=0.9) +
  geom_line(aes(y = (CRWSST), colour = "CRW_SST"), linetype=1, size=0.3, alpha=0.9) +
   geom_line(aes(y = (SST_CCI), colour = "SST_CCI"), size=0.4, alpha=0.9) +
  
  scale_x_date(date_breaks = "3 months", date_labels = "%y-%m", 
               limits = c(as.Date("2014-07-15"), as.Date("2016-07-10")),
               name = "") +
  scale_y_continuous(limits = c(26,31.1), expand = c(0.01, 0), 
                   name=(expression("Temperature"~(degree*C))),
                   breaks = seq(26,31, by=1))+
  annotate("rect", xmin= as.Date("2014-08-17"), xmax = as.Date("2014-08-24"), 
                 ymin = 26, ymax = 31, fill="blue", alpha = .3) +
  annotate("rect", xmin= as.Date("2015-08-16"), xmax = as.Date("2015-08-22"), 
                 ymin = 26, ymax = 31, fill="blue", alpha = .3) + 
  annotate("rect", xmin= as.Date("2016-04-16"), xmax = as.Date("2016-04-21"), 
                 ymin = 26, ymax = 31, fill="blue", alpha = .3)
  
DHWPlot<-ggplot(DHW.Wide, aes(x=Date)) + ggthe_bw+
   theme(legend.box.background = element_rect(),
        legend.position = c(0.2, 0.5)) +
    scale_colour_manual(values = c("red", "gray", "black", "purple"), name="Source")
DHWPlot<-DHWPlot +
    geom_hline(yintercept=8, linetype="dashed", color="red", size=0.5)+
    geom_hline(yintercept=4, linetype="dashed", color="orange", size=0.5)+
    geom_line(aes(y = (OI_DHW), colour = "OI SST"), linetype=1, size=0.3, alpha=0.9) +      
    geom_line(aes(y = (InSitu_DHW), colour = "In situ"), size=0.4, alpha=0.9) +
    geom_line(aes(y = (CRW_DHW), colour = "CRW"), size=0.4, alpha=0.9) +  
    geom_line(aes(y = (CCI_DHW), colour = "CCI"), size=0.4, alpha=0.9) +  
    
    scale_x_date(date_breaks = "3 months", date_labels = "%m/%Y", 
               limits = c(as.Date("2014-07-15"), as.Date("2016-07-10")),
               name = "Month-Year") +
    scale_y_continuous(limits = c(0,10), expand = c(0.01, 0), 
                   name=(expression("DHW"~(degree*C~"-weeks"))),
                   breaks = seq(0,10, by=2))+
    annotate("rect", xmin= as.Date("2014-08-17"), xmax = as.Date("2014-08-24"), 
                 ymin = 0, ymax = 10, fill="blue", alpha = .3) +
    annotate("rect", xmin= as.Date("2015-08-16"), xmax = as.Date("2015-08-22"), 
                 ymin = 0, ymax = 10, fill="blue", alpha = .3) + 
    annotate("rect", xmin= as.Date("2016-04-16"), xmax = as.Date("2016-04-21"), 
                 ymin = 0, ymax = 10, fill="blue", alpha = .3)

Figure2AllSources<-grid.arrange(TemperaturePlotAll, DHWPlot, nrow=2, heights=c(6/10, 4/10))

#ggsave(file="C.Outputs/Figure_1b.svg", plot=Figure2a_Final, width=6.5, height=4.5)

Figure 2 (Supplementary): Accumulated heat stress and thermal regime at Uva Island, Gulf of Chiriquí. a. Thermal regime before and during the 2015-16 El Niño using daily OISST (black) and in situ (gray) data collected at Uva Island reef (~3m depth). Upper plot and right-hand axis show mean daily temperature (solid lines) and bleaching threshold (dashed lines) for each data set. Lower plot and left-hand axis show DHW calculated from each data source. The blue vertical lines represent the periods when the coral samples were collected.
# 5. Potential bias based on data source (Supplementary Information)

5.1 SST with data from 2014-2016

  Temperature.Wide$Year<-year(Temperature.Wide$Date)
  Temperature.Wide_2014_2016<-filter(Temperature.Wide, Year>2013)
  Temperature.Wide_2014_2016<-filter(Temperature.Wide_2014_2016, Year<2017)
  Temperature.Wide_2014_2016<-Temperature.Wide_2014_2016 %>% drop_na(Temperature)
  summary(Temperature.Wide_2014_2016$OISST)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   26.82   28.73   29.25   29.20   29.77   30.84
  summary(Temperature.Wide_2014_2016$CRWSST )
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   27.28   28.39   28.91   28.88   29.36   30.44
  summary(Temperature.Wide_2014_2016$Temperature)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   26.70   28.63   29.17   29.15   29.81   30.78
  summary(Temperature.Wide_2014_2016$SST_CCI)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   27.06   28.78   29.12   29.18   29.68   30.81
  cov_OISST<-cov(Temperature.Wide_2014_2016$Temperature, Temperature.Wide_2014_2016$OISST)
  cov_CCI<-cov(Temperature.Wide_2014_2016$Temperature, Temperature.Wide_2014_2016$SST_CCI )
  cov_CRW<-cov(Temperature.Wide_2014_2016$Temperature, Temperature.Wide_2014_2016$CRWSST)
  
  cor_OISST<-cor(Temperature.Wide_2014_2016$Temperature, Temperature.Wide_2014_2016$OISST)
  cor_CCI<-cor(Temperature.Wide_2014_2016$Temperature, Temperature.Wide_2014_2016$SST_CCI )
  cor_CRW<-cor(Temperature.Wide_2014_2016$Temperature, Temperature.Wide_2014_2016$CRWSST)

Figure S6 - OI, CRW

## 
## Call:
## lm(formula = OISST ~ Temperature, data = Temperature.Wide_2014_2016)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3.10534 -0.24675  0.03016  0.28238  1.48478 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.90978    0.53350   12.95   <2e-16 ***
## Temperature  0.76482    0.01829   41.80   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4323 on 866 degrees of freedom
## Multiple R-squared:  0.6687, Adjusted R-squared:  0.6683 
## F-statistic:  1748 on 1 and 866 DF,  p-value: < 2.2e-16
## 
## Call:
## lm(formula = CRWSST ~ Temperature, data = Temperature.Wide_2014_2016)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.06653 -0.31235 -0.06539  0.27989  1.58472 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 10.42236    0.53621   19.44   <2e-16 ***
## Temperature  0.63324    0.01839   34.44   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4345 on 866 degrees of freedom
## Multiple R-squared:  0.578,  Adjusted R-squared:  0.5775 
## F-statistic:  1186 on 1 and 866 DF,  p-value: < 2.2e-16

Supplementary Figure S6: Temperature comparison among different data sources. a. Temperature for the August 2014 - August 2016 period. b. Correlation between OISST (v2) data and in situ (3m) data. c. Correlation between CRW (CoralTemp v3.1) and in situ data. Panels I contain SST values under the bleaching threshold for both data sets being compared (there is no accumulation of DHW). Panels II contain SST values above the bleaching threshold for both data sets being compared (both data sets accumulate of DHW). Panels III contain the cases in which satellite values were above the bleaching threshold, both not the in situ values (potential DHW overestimation). Panels IV contain the cases in which in situ values were above the bleaching threshold, both not the satellite values (potential DHW underestimation).

Figure S6 - OI, CRW, CCI - Final MS

## 
## Call:
## lm(formula = OISST ~ Temperature, data = Temperature.Wide_2014_2016)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3.10534 -0.24675  0.03016  0.28238  1.48478 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.90978    0.53350   12.95   <2e-16 ***
## Temperature  0.76482    0.01829   41.80   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4323 on 866 degrees of freedom
## Multiple R-squared:  0.6687, Adjusted R-squared:  0.6683 
## F-statistic:  1748 on 1 and 866 DF,  p-value: < 2.2e-16
## 
## Call:
## lm(formula = CRWSST ~ Temperature, data = Temperature.Wide_2014_2016)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.06653 -0.31235 -0.06539  0.27989  1.58472 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 10.42236    0.53621   19.44   <2e-16 ***
## Temperature  0.63324    0.01839   34.44   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4345 on 866 degrees of freedom
## Multiple R-squared:  0.578,  Adjusted R-squared:  0.5775 
## F-statistic:  1186 on 1 and 866 DF,  p-value: < 2.2e-16
## 
## Call:
## lm(formula = SST_CCI ~ Temperature, data = Temperature.Wide_2014_2016)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.53788 -0.21328 -0.00374  0.22694  1.41952 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 11.43513    0.46270   24.71   <2e-16 ***
## Temperature  0.60866    0.01587   38.36   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3749 on 866 degrees of freedom
## Multiple R-squared:  0.6295, Adjusted R-squared:  0.6291 
## F-statistic:  1471 on 1 and 866 DF,  p-value: < 2.2e-16
## 
## Call:
## lm(formula = OI_HS ~ IS_HS, data = Temperature.Wide_2014_2016)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.78617 -0.19605  0.00501  0.20607  0.95227 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.18192    0.03148    5.78  1.5e-08 ***
## IS_HS        0.62786    0.03709   16.93  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2969 on 401 degrees of freedom
##   (465 observations deleted due to missingness)
## Multiple R-squared:  0.4168, Adjusted R-squared:  0.4153 
## F-statistic: 286.6 on 1 and 401 DF,  p-value: < 2.2e-16
## 
## Call:
## lm(formula = CRW_HS ~ IS_HS, data = Temperature.Wide_2014_2016)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.63902 -0.24560 -0.01406  0.24626  0.79704 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.25274    0.03401    7.43 7.84e-13 ***
## IS_HS        0.41891    0.03907   10.72  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3083 on 363 degrees of freedom
##   (503 observations deleted due to missingness)
## Multiple R-squared:  0.2406, Adjusted R-squared:  0.2385 
## F-statistic:   115 on 1 and 363 DF,  p-value: < 2.2e-16
## 
## Call:
## lm(formula = CCI_HS ~ IS_HS, data = Temperature.Wide_2014_2016)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.61703 -0.17557 -0.02269  0.16680  0.93975 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.17896    0.03356   5.333 1.77e-07 ***
## IS_HS        0.43455    0.03679  11.811  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2539 on 337 degrees of freedom
##   (529 observations deleted due to missingness)
## Multiple R-squared:  0.2928, Adjusted R-squared:  0.2907 
## F-statistic: 139.5 on 1 and 337 DF,  p-value: < 2.2e-16

Supplementary Figure S6: Temperature comparison among different data sources. a. Temperature for the August 2014 - August 2016 period. b. Correlation between OISST (v2) data and in situ (3m) data. c. Correlation between CRW (CoralTemp v3.1) and in situ data. d. Correlation between CCI (ESA v2.1) and in situ data. Panels I contain SST values under the bleaching threshold for both data sets being compared (there is no accumulation of DHW). Panels II contain SST values above the bleaching threshold for both data sets being compared (both data sets accumulate of DHW). Panels III contain the cases in which satellite values were above the bleaching threshold, both not the in situ values (potential DHW overestimation). Panels IV contain the cases in which in situ values were above the bleaching threshold, both not the satellite values (potential DHW underestimation). Panels e., f. and g. Show the correlation between the Hotspots (HS) for each satellite source compared to in situ data.

5.2 SST with all data from any year with the 3 sources

Figure OI, CRW

  Temperature.Wide.All<-Temperature.Wide[complete.cases(Temperature.Wide), ]
  cov_OISST_all<-cov(Temperature.Wide.All$Temperature, Temperature.Wide.All$OISST)
  cov_CCI_all<-cov(Temperature.Wide.All$Temperature, Temperature.Wide.All$SST_CCI )
  cov_CRW_all<-cov(Temperature.Wide.All$Temperature, Temperature.Wide.All$CRWSST)
  
  cor_OISST_all<-cor(Temperature.Wide.All$Temperature, Temperature.Wide.All$OISST)
  cor_CCI_all<-cor(Temperature.Wide.All$Temperature, Temperature.Wide.All$SST_CCI )
  cor_CRW_all<-cor(Temperature.Wide.All$Temperature, Temperature.Wide.All$CRWSST)
## 
## Call:
## lm(formula = OISST ~ Temperature, data = Temperature.Wide.All)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3.07279 -0.31868  0.04338  0.34622  1.99901 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 5.826754   0.256110   22.75   <2e-16 ***
## Temperature 0.799726   0.008955   89.30   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5533 on 5813 degrees of freedom
## Multiple R-squared:  0.5784, Adjusted R-squared:  0.5783 
## F-statistic:  7975 on 1 and 5813 DF,  p-value: < 2.2e-16
## 
## Call:
## lm(formula = CRWSST ~ Temperature, data = Temperature.Wide.All)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.57787 -0.26038 -0.00314  0.24734  1.97499 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 7.916719   0.188561   41.98   <2e-16 ***
## Temperature 0.719924   0.006593  109.19   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4073 on 5813 degrees of freedom
## Multiple R-squared:  0.6722, Adjusted R-squared:  0.6722 
## F-statistic: 1.192e+04 on 1 and 5813 DF,  p-value: < 2.2e-16
## 
## Call:
## lm(formula = SST_CCI ~ Temperature, data = Temperature.Wide.All)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.53491 -0.24624 -0.01689  0.22996  1.93370 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 9.781716   0.181576   53.87   <2e-16 ***
## Temperature 0.666116   0.006349  104.92   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3922 on 5813 degrees of freedom
## Multiple R-squared:  0.6544, Adjusted R-squared:  0.6544 
## F-statistic: 1.101e+04 on 1 and 5813 DF,  p-value: < 2.2e-16

Figure OI, CRW, CCI

## 
## Call:
## lm(formula = OISST ~ Temperature, data = Temperature.Wide.All)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3.07279 -0.31868  0.04338  0.34622  1.99901 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 5.826754   0.256110   22.75   <2e-16 ***
## Temperature 0.799726   0.008955   89.30   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5533 on 5813 degrees of freedom
## Multiple R-squared:  0.5784, Adjusted R-squared:  0.5783 
## F-statistic:  7975 on 1 and 5813 DF,  p-value: < 2.2e-16
## 
## Call:
## lm(formula = CRWSST ~ Temperature, data = Temperature.Wide.All)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.57787 -0.26038 -0.00314  0.24734  1.97499 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 7.916719   0.188561   41.98   <2e-16 ***
## Temperature 0.719924   0.006593  109.19   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4073 on 5813 degrees of freedom
## Multiple R-squared:  0.6722, Adjusted R-squared:  0.6722 
## F-statistic: 1.192e+04 on 1 and 5813 DF,  p-value: < 2.2e-16
## 
## Call:
## lm(formula = SST_CCI ~ Temperature, data = Temperature.Wide.All)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.53491 -0.24624 -0.01689  0.22996  1.93370 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 9.781716   0.181576   53.87   <2e-16 ***
## Temperature 0.666116   0.006349  104.92   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3922 on 5813 degrees of freedom
## Multiple R-squared:  0.6544, Adjusted R-squared:  0.6544 
## F-statistic: 1.101e+04 on 1 and 5813 DF,  p-value: < 2.2e-16
## 
## Call:
## lm(formula = OI_HS ~ IS_HS, data = Temperature.Wide.All)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.83251 -0.22729 -0.02984  0.21192  1.21630 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.25431    0.01687   15.07   <2e-16 ***
## IS_HS        0.54762    0.02412   22.70   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3231 on 1097 degrees of freedom
##   (4716 observations deleted due to missingness)
## Multiple R-squared:  0.3196, Adjusted R-squared:  0.319 
## F-statistic: 515.3 on 1 and 1097 DF,  p-value: < 2.2e-16
## 
## Call:
## lm(formula = CRW_HS ~ IS_HS, data = Temperature.Wide.All)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.73740 -0.19661 -0.01626  0.17122  1.05972 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.20148    0.01416   14.22   <2e-16 ***
## IS_HS        0.52742    0.02065   25.54   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2836 on 1161 degrees of freedom
##   (4652 observations deleted due to missingness)
## Multiple R-squared:  0.3598, Adjusted R-squared:  0.3592 
## F-statistic: 652.4 on 1 and 1161 DF,  p-value: < 2.2e-16
## 
## Call:
## lm(formula = CCI_HS ~ IS_HS, data = Temperature.Wide.All)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.67164 -0.18725 -0.04442  0.14948  1.26001 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.17748    0.01525   11.64   <2e-16 ***
## IS_HS        0.47954    0.02096   22.88   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.268 on 1042 degrees of freedom
##   (4771 observations deleted due to missingness)
## Multiple R-squared:  0.3343, Adjusted R-squared:  0.3337 
## F-statistic: 523.3 on 1 and 1042 DF,  p-value: < 2.2e-16

5.3 DHW with data from 2014-2016

# Organize the data source
Wide.data.DHW<-dplyr::filter(Long.data.DHW, Year>2013)
Wide.data.DHW<-dplyr::filter(Wide.data.DHW, Year<2017)
Wide.data.DHW<-Wide.data.DHW %>% 
              dplyr::select(Source, Date, DHW)

Wide.data.DHW2<-reshape(Wide.data.DHW, 
                        idvar = "Date", timevar = "Source", direction = "wide")
Wide.data.DHW2<-Wide.data.DHW2[complete.cases (Wide.data.DHW2),]
  cov_OIDHW<-cov(Wide.data.DHW2$DHW.In_situ, Wide.data.DHW2$DHW.OI_SST)
  cov_CCIDHW<-cov(Wide.data.DHW2$DHW.In_situ, Wide.data.DHW2$DHW.CCI_SST)
  cov_CRWDHW<-cov(Wide.data.DHW2$DHW.In_situ, Wide.data.DHW2$DHW.CRW_SST)
  
  cor_OIDHW<-cor(Wide.data.DHW2$DHW.In_situ, Wide.data.DHW2$DHW.OI_SST)
  cor_CCIDHW<-cor(Wide.data.DHW2$DHW.In_situ, Wide.data.DHW2$DHW.CCI_SST)
  cor_CRWDHW<-cor(Wide.data.DHW2$DHW.In_situ, Wide.data.DHW2$DHW.CRW_SST)


#head(Wide.data.DHW2)
#tail(Wide.data.DHW2)

Figure S7

## 
## Call:
## lm(formula = DHW.OI_SST ~ DHW.In_situ, data = Wide.data.DHW2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.1786 -0.2485  0.3756  0.3756  1.6963 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.375601   0.034351  -10.93   <2e-16 ***
## DHW.In_situ  0.915716   0.009934   92.18   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7367 on 783 degrees of freedom
## Multiple R-squared:  0.9156, Adjusted R-squared:  0.9155 
## F-statistic:  8496 on 1 and 783 DF,  p-value: < 2.2e-16
## 
## Call:
## lm(formula = DHW.CRW_SST ~ DHW.In_situ, data = Wide.data.DHW2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.83582 -0.23727 -0.08136  0.38007  2.98418 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.08136    0.05460    1.49    0.137    
## DHW.In_situ  0.37167    0.01579   23.54   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.171 on 783 degrees of freedom
## Multiple R-squared:  0.4144, Adjusted R-squared:  0.4137 
## F-statistic: 554.1 on 1 and 783 DF,  p-value: < 2.2e-16
## 
## Call:
## lm(formula = DHW.CCI_SST ~ DHW.In_situ, data = Wide.data.DHW2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.67425 -0.21116 -0.07320  0.00814  1.15500 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.073202   0.019004   3.852 0.000127 ***
## DHW.In_situ 0.123145   0.005496  22.406  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4076 on 783 degrees of freedom
## Multiple R-squared:  0.3907, Adjusted R-squared:  0.3899 
## F-statistic:   502 on 1 and 783 DF,  p-value: < 2.2e-16

Supplementary Figure S7: Degree heating weeks (DHW) index comparison among different data sources. a. DHW calculated for November 2014 - August 2016. b. Correlation between OISST derived in situ derived DHW indices. c. Correlation between CRW (CoralTemp v3.1) derived and in situ derived DHW indices. d. Correlation between CIISST derived in situ derived DHW indices.

5.4 DHW with all DHW data available for all 3 sources

# Organize the data source

  Wide.data.DHW<-Long.data.DHW %>% 
                dplyr::select(Source, Date, DHW)
  
  Wide.data.DHW3<-reshape(Wide.data.DHW, 
                          idvar = "Date", timevar = "Source", direction = "wide")
  
  # Wide.data.DHW2<-Wide.data.DHW2[complete.cases (Wide.data.DHW2),]
  #head(Wide.data.DHW2)
  #tail(Wide.data.DHW2)
# a. Estimated DHW

    DHWSource2<-ggplot(Wide.data.DHW3, aes(x=Date)) + 
       ggthe_bw + ggtitle("a.") +
       scale_x_date(date_breaks = "24 months", date_labels = "%Y", 
                    limits = c(as.Date("1997-01-01"), as.Date("2016-12-31")),
                    expand = c(0.01, 0), 
                    name = "Month-Year") +
      scale_y_continuous(limits = c(0,11), expand = c(0.01, 0), 
                         name=(expression("DHW"~(degree*C~-"weeks"))),
                         breaks = seq(0,10, by=2))+
      theme(legend.box.background = element_rect(),
            legend.position = c(0.25, 0.7))+
      
      geom_line(aes(y = (DHW.OI_SST), colour = "OISST"), linetype=1, size=0.3, alpha=0.9) +      
      geom_line(aes(y = (DHW.In_situ), colour = "In situ"), size=0.4, alpha=0.9) +
      geom_line(aes(y = (DHW.CCI_SST), colour = "CCI"), size=0.4, alpha=0.9) +
      geom_line(aes(y = (DHW.CRW_SST), colour = "CRW"), linetype=1, size=0.3, alpha=0.9) +
      scale_colour_manual(values = c("red", "blue4","grey","black"), name=" Data source")
    #DHWSource2
    #ggsave(file="Outputs/Sup_DHW_Source2.svg", plot=DHWSource2, width=8, height=3.5)


# b. OI versus In situ 

    LM_OI_DWW2<-lm(DHW.OI_SST~DHW.In_situ, data = Wide.data.DHW3)
    summary(LM_OI_DWW2)
## 
## Call:
## lm(formula = DHW.OI_SST ~ DHW.In_situ, data = Wide.data.DHW3)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.44465 -0.06661 -0.06661 -0.06661  2.06520 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.066607   0.005964   11.17   <2e-16 ***
## DHW.In_situ 0.852559   0.003937  216.57   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4174 on 5398 degrees of freedom
##   (8093 observations deleted due to missingness)
## Multiple R-squared:  0.8968, Adjusted R-squared:  0.8968 
## F-statistic: 4.69e+04 on 1 and 5398 DF,  p-value: < 2.2e-16
    GCITYvsOI_DHW2<-ggplot(Wide.data.DHW3, aes(x=DHW.In_situ, y=DHW.OI_SST)) + 
      ggthe_bw + ggtitle("b.")+
      
      geom_abline(slope = 1, intercept = 0, colour="grey")+
      geom_point(alpha=0.5, size=1)+ 
      geom_smooth(method = lm, colour="black")+
      
      scale_y_continuous(limits = c(0,9), expand = c(0.01, 0), 
                         name=(expression("OISST DHW"~(degree*C~-weeks))),
                         breaks = seq(0,9, by=1))+
      scale_x_continuous(limits = c(0,9), expand = c(0.01, 0), 
                         name=(expression("In situ DHW"~(degree*C~-weeks))),
                         breaks = seq(0,9, by=1))+
     
      annotate("text", x = 4.2, y = 8.2, 
               label = "OISST_DHW = 0.8 (In situ_DHW) + 0.1 \n R2: 0.9",
               size=3)
    #GCITYvsOI_DHW2
    # ggsave(file="Outputs/DHW_IS_OI_2.svg", plot=GCITYvsOI_DHW2, width=4, height=3.5)


# c. CRW versus In situ 

    LM_CRW_DHW2<-lm(DHW.CRW_SST~DHW.In_situ, data = Wide.data.DHW3)
    summary(LM_CRW_DHW2)
## 
## Call:
## lm(formula = DHW.CRW_SST ~ DHW.In_situ, data = Wide.data.DHW3)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.4617 -0.0452 -0.0452 -0.0452  3.5324 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.045180   0.008011    5.64 1.79e-08 ***
## DHW.In_situ 0.461016   0.005288   87.19  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5607 on 5398 degrees of freedom
##   (8093 observations deleted due to missingness)
## Multiple R-squared:  0.5848, Adjusted R-squared:  0.5847 
## F-statistic:  7602 on 1 and 5398 DF,  p-value: < 2.2e-16
    GCITYvsCRW_DHW2<-ggplot(Wide.data.DHW3, aes(x=DHW.In_situ, y=DHW.CRW_SST)) + 
      ggthe_bw + ggtitle("c.")+
      
      geom_abline(slope = 1, intercept = 0, colour="grey")+ 
      geom_point(alpha=0.5, size=1)+ 
      geom_smooth(method = lm, colour="blue4")+
      
      scale_y_continuous(limits = c(0,9), expand = c(0.01, 0), 
                         name=(expression("CRW DHW"~(degree*C~-weeks))),
                         breaks = seq(0,10, by=1))+
      scale_x_continuous(limits = c(0,9), expand = c(0.01, 0), 
                         name=(expression("In situ DHW "~(degree*C~-weeks))),
                         breaks = seq(0,10, by=1))+
      annotate("text", x = 4.2, y = 8.2, 
               label = "CRW_DHW = 0.5 (In situ_DHW) \n R2: 0.6",
               size=3)
    #GCITYvsCRW_DHW2
    #ggsave(file="Outputs/DHW_IS_CRW_2.svg", plot=GCITYvsCRW_DHW2, width=4, height=3.5)


# d. CCI versus In situ 

    LM_CCI_DHW2<-lm(DHW.CCI_SST~DHW.In_situ, data = Wide.data.DHW3)
    summary(LM_CCI_DHW2)
## 
## Call:
## lm(formula = DHW.CCI_SST ~ DHW.In_situ, data = Wide.data.DHW3)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.45379 -0.09246 -0.09246 -0.09246  2.31778 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.092464   0.006510   14.20   <2e-16 ***
## DHW.In_situ 0.227753   0.004297   53.01   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4556 on 5398 degrees of freedom
##   (8093 observations deleted due to missingness)
## Multiple R-squared:  0.3423, Adjusted R-squared:  0.3422 
## F-statistic:  2810 on 1 and 5398 DF,  p-value: < 2.2e-16
    GCITYvsCCI_DHW2<-ggplot(Wide.data.DHW3, aes(x=DHW.In_situ, y=DHW.CCI_SST)) + 
      ggthe_bw + ggtitle("d.")+
      
      geom_abline(slope = 1, intercept = 0, colour="grey")+ 
      geom_point(alpha=0.5, size=1)+ 
      geom_smooth(method = lm, colour="red")+
      
      scale_y_continuous(limits = c(0,9), expand = c(0.01, 0), 
                         name=(expression("CRW DHW"~(degree*C~-weeks))),
                         breaks = seq(0,10, by=1))+
      scale_x_continuous(limits = c(0,9), expand = c(0.01, 0), 
                         name=(expression("In situ DHW "~(degree*C~-weeks))),
                         breaks = seq(0,10, by=1))+
      annotate("text", x = 4.2, y = 8.2, 
               label = "CCI_DHW = 0.2 (In situ_DHW) + 0.1 \n R2: 0.3",
               size=3)
    #GCITYvsCCI_DHW2
    #ggsave(file="Outputs/DHW_IS_CRW_2.svg", plot=GCITYvsCRW_DHW2, width=4, height=3.5)


# d. OI differences 

    Wide.data.DHW3$OI_Difference<-
      (Wide.data.DHW3$DHW.OI_SST-Wide.data.DHW3$DHW.In_situ)
    
    OI_DHW_Residuals2<-ggplot(Wide.data.DHW3, aes(x=DHW.In_situ, y=OI_Difference)) + 
      ggthe_bw+ ggtitle("d.") +
      
      geom_abline(slope = 0, intercept = 0, colour="grey") +
      geom_point(alpha=0.5, size=1)+ 
      geom_smooth(method = lm, colour="black")+
      
      scale_y_continuous(limits = c(-3,3), expand = c(0, 0), 
                         name=(expression("OI difference"~(degree*C~-weeks))),
                         breaks = seq(-3,3, by=0.5))+
      scale_x_continuous(limits = c(0,9), expand = c(0.01, 0), 
                         name=(expression("In situ DHW"~(degree*C~-weeks))),
                         breaks = seq(0,10, by=1))
    #OI_DHW_Residuals2
    #ggsave(file="Outputs/DHW_OI_Differences.svg", plot=OI_DHW_Residuals2, width=4, height=3.5)

    
# e. OI differences 

    Wide.data.DHW3$CRW_Difference<-
      (Wide.data.DHW3$DHW.CRW_SST-Wide.data.DHW3$DHW.In_situ)
    
    CRW_DHW_Residuals2<-ggplot(Wide.data.DHW3, aes(x=DHW.In_situ, y=CRW_Difference)) + 
      ggthe_bw + ggtitle("e.") +
      
      geom_abline(slope = 0, intercept = 0, colour="grey")+ 
      geom_point(alpha=0.5, size=1)+ 
      geom_smooth(method = lm, colour="blue4")+
      
      scale_y_continuous(limits = c(-3,3), expand = c(0, 0), 
                         name=(expression("CRW difference"~(degree*C))),
                         breaks = seq(-3,3, by=0.5))+
      scale_x_continuous(limits = c(0,9), expand = c(0.01, 0), 
                         name=(expression("In situ"~(degree*C))),
                         breaks = seq(0,10, by=1))
    #CRW_DHW_Residuals2
    #ggsave(file="Outputs/DHW_CRW_Differences.svg", plot=CRW_DHW_Residuals2, width=4, height=3.5)


Sources4<-grid.arrange(
  DHWSource2, arrangeGrob(GCITYvsOI_DHW2,GCITYvsCRW_DHW2, GCITYvsCCI_DHW2, ncol=3), heights=c(5/10, 5/10), ncol = 1)

#ggsave(file="Outputs/DHW_Sourcess2.svg", plot=Sources4, width=7, height=8)

5.5 CCI vs OI DHW all ENSOS

Figure S8

FigureS6a<-ggplot(Nino.OI_data, aes(x=Days)) + ggtitle("a.")+ ggthe_bw+
  theme(legend.box.background = element_rect(),
        legend.position = c(0.25, 0.75),
         panel.grid= element_blank()) +
  geom_hline(yintercept=8, linetype="dashed", color="red", size=1)+
  geom_hline(yintercept=4, linetype="dashed", color="orange", size=1)+
  geom_line(aes(y = DHW, colour = factor(Event)), size=1.0) + 
  scale_x_continuous(limits = c(0,701), expand = c(0, 0), 
                     name=("Month"),
                     breaks = seq(0,701, by= 90))+
  scale_y_continuous(limits = c(0,11), 
                     name=(expression("DHW"~(degree*C ~"-weeks"))), 
                     expand = c(0.01, 0.01),
                     breaks = seq(0,10, by=2)) +  
  scale_colour_manual(values = c("magenta","cyan", "black"), name="ENSO")
#FigureS5a

FigureS6b<-ggplot(Nino.CCI_data, aes(x=Days)) + ggtitle("b.")+ ggthe_bw+
  theme(legend.box.background = element_rect(),
        legend.position = c(0.25, 0.75),
         panel.grid= element_blank()) +
  geom_hline(yintercept=8, linetype="dashed", color="red", size=1)+
  geom_hline(yintercept=4, linetype="dashed", color="orange", size=1)+
  geom_line(aes(y = DHW, colour = factor(Event)), size=1.0) + 
  scale_x_continuous(limits = c(0,701), expand = c(0, 0), 
                     name=("Month"),
                     breaks = seq(0,701, by= 90))+
  scale_y_continuous(limits = c(0,11), 
                     name=(expression("DHW"~(degree*C ~"-weeks"))), expand = c(0.01, 0.01),
                     breaks = seq(0,10, by=2)) +  
  scale_colour_manual(values = c("magenta","cyan", "black"), name="ENSO")
#FigureS5b

Sources6<-grid.arrange(FigureS6a,FigureS6b, ncol=2)

#ggsave(file="C.Outputs/Figure_S5.svg", plot=Sources6, width=8, height=4.5)

6. Other ENSO characteristics (Supplementary Information)

6.1 HotSpots magnitude

Other way to compare ENSOs is to look at the magnitude and number of the Hot Spots

Suplemantary Figure S9 - HotSpots

# Supplemantary Figure S9
  
  SupplumentFigure3a<-ggplot(Nino.OI_data, aes(x=Days)) + ggthe_bw +
      theme(legend.box.background = element_rect(),
            legend.position = "none") +
      
      geom_hline(yintercept=1, linetype="dashed", color="black")+
      geom_hline(yintercept=2, linetype="dashed", color="grey")+
      scale_x_continuous(limits = c(0,701), expand = c(0, 0), 
                         name=("Month of the year"),
                         breaks = seq(0,701, by= 90))+
      scale_y_continuous(limits = c(0,3), 
                       name=(expression("HotSpot"~(degree*C))), expand = c(0.02, 0.03),
                         breaks = seq(0,3, by=0.5))
  
  SupplumentFigure3a <-SupplumentFigure3a +
      #geom_point(aes(y = HotSpot, colour = factor(Event)), size=0.02) + 
      geom_line(aes(y = HotSpot, colour = factor(Event))) + 
     
      scale_colour_manual(values = c("magenta","cyan", "black"), name="ENSO")+
      #geom_point(aes(y = Clima_Norm, size=1)) +
      facet_grid(Event~.)
  SupplumentFigure3a

  #ggsave(file="C.Outputs/Figure_S3.svg", plot=SupplumentFigure3a, width=6, height=5)

Supplementary Figure S9: HotSpots (SST - MMM) at Uva Island reef during the two-year period of each El Niño calculated using NOAA OISSTv2 data. Number of HotSpots (n) have been summarized to show how many times they exceeded 2ºC (HotSpots > 2ºC, top values in each panel), were higher than 1ºC but lower or equal to 2ºC (1ºC < HotSpots ≤ 2ºC, middle values), or were lower or equal to 1ºC (0ºC < HotSpots ≤ 1ºC, bottom values). Note that only HotSpots > 1ºC are included in the estimation of DHW. Summation of all HotSpots > 1ºC during each 2 year period were 131.2 during 1982-83, 123.5 during 1997-98, and 108.7 during 2015-16.

6.2 Number of HS higher than 0, 1, 2, and 3C

Number of HS>1 (used in DHW calculation)

    # Number of HotSpots; HS > 1
        HS_N_1<-as.data.frame(Nino.OI_data %>%
              group_by(Event) %>%
              summarize(N = n(),
                        Cont=sum (HotSpot>=1)))
        HS_N_1

Number of 0<HS<1 (not used in DHW calculation)

# Number of HotSpots; 0<HS<1
        HS_N_0_1<-as.data.frame(Nino.OI_data %>%
              group_by(Event) %>%
              summarize(N = n(),
                        Cont=sum (HotSpot>0 & HotSpot<1)))
        HS_N_0_1

Number of 1<HS<2 (used in DHW calculation, lower stress)

 # Number of HotSpots; 1<HS<2
        HS_N_1_2<-as.data.frame(Nino.OI_data %>%
              group_by(Event) %>%
              summarize(N = n(),
                        Cont=sum (HotSpot>=1 & HotSpot<2)))
        HS_N_1_2

Number of 2<HS<3 (used in DHW calculation, higher stress)

# Number of HotSpots; 2<HS<3
        HS_N_2_3<-as.data.frame(Nino.OI_data %>%
              group_by(Event) %>%
              summarize(N = n(),
                        Cont=sum (HotSpot>=2 & HotSpot<3)))
        HS_N_2_3

6.3 Area under HS curve

Area under the HS curve (all HS added)

# Area under of (HS>0)
  HS_Aggregated<-as.data.frame(Nino.OI_data %>%
        group_by(Event) %>%
        summarize(Sum = sum (HotSpot)))
  HS_Aggregated

Area under the HS curve HS>1 added (used in DHW)

# Area under of (HS>1)
  HS1_Aggregated<-as.data.frame(Nino.OI_data %>%
        group_by(Event) %>%
        summarize(Sum = sum (HotSpot[HotSpot>1])))
  HS1_Aggregated

Area under the HS curve only HS>2 added

# Area under (HS>2)
  HS2_Aggregated<-as.data.frame(Nino.OI_data %>%
        group_by(Event) %>%
        summarize(Sum = sum (HotSpot[HotSpot>2])))
  HS2_Aggregated

7. Future DHW

#SSP5-8.5_SST (provided by Ruben van Hooidonk)
    Uva.Future<-read.csv("B.Temperature_data/Uva_DHW_SSP5-8.5_All.csv",header = T)
head(Uva.Future)
Uva.Future$Difference<-Uva.Future$Max.DHW-Uva.Future$Max.DHW.MMM.1

Figure 5b: Maximum DHW projected for the period of 2018-2100 at Uva Island reef

DHW_future<-ggplot(Uva.Future) + 
      geom_line(aes(x=Year, y=Max.DHW.MMM.1, color="SSP 5-8.8"), 
                linetype=1, size=0.3, alpha=1)+
  geom_line(aes(x=Year, y=Max.DHW.MMM.2.25, color="SSP 5-8.8 + 1.25 C tolerance"), 
                linetype=1, size=0.3, alpha=1)+
      geom_line(aes(x=Year, y=Max.DHW.MMM.2.5, color="SSP 5-8.8 + 1.50 C tolerance"), 
                linetype=1, size=0.3, alpha=1)+
  scale_color_manual(name = "Scenario", values = c("SSP 5-8.8" = "red", 
                                                   "SSP 5-8.8 + 1.25 C tolerance" = "purple", 
                                                   "SSP 5-8.8 + 1.50 C tolerance" = "blue"))+
  theme_bw()+
  theme(legend.box.background = element_rect(),
        legend.position = "top",
        panel.grid= element_blank(), 
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  
   geom_hline(yintercept=4.00, linetype="dashed", 
              color="orange", size=0.5)+
   geom_hline(yintercept=8.00, linetype="dashed",
              color="red", size=0.5)+
  
  # geom_hline(yintercept=10.00, linetype="dashed", 
  #            color="black", size=0.5)+
  # geom_hline(yintercept=20.00, linetype="dashed",
  #            color="black", size=1)+
  # geom_hline(yintercept=30.00, linetype="dashed",
  #            color="black", size=1.5)+
  # geom_hline(yintercept=40.00, linetype="dashed",
  #            color="black", size=2)+
  
  scale_x_continuous(limits = c(2018, 2100),
                     name=(""),
                     breaks = seq(2020, 2100, by= 5),
                    expand = c(0.01, 0.01))+
  scale_y_continuous(limits = c(0,42), 
                     name=(expression(
                       "Expected annual maximum DHW"~(
                         degree*C ~"-weeks")~"under SSP-8.5")), 
                     expand = c(0.01, 0.01),
                     breaks = seq(0,40, by=4))#+
 # annotate("text", x = c(2083, 2030, 2030, 2030),
 #                  y = c(11, 21, 31, 41), 
 #          label = c("Max DHW recorded (1998)",
 #                    "2x max DHW recorded",
 #                    "3x max DHW recorded", 
 #                    "4x max DHW recorded"))
  
DHW_future

#ggsave(file="C.Outputs/Figure_SSP.svg", plot=DHW_future, width=5.5, height=5.0)

Figure 5b: Maximum DHW projected for the period of 2018-2100 at Uva Island reef using Coupled Model Intercomparison Projects climate models (CMIP6) under the Shared Socioeconomic Pathways (SSP) 5-8.5 scenario and three different bleaching thresholds (MMM+1C, MMM+2.25C and MMM+2.5C).