ESR_ACLIM.rmd

---
title: "ACLIM2 ESR indices"
author: K. Holsman
always_allow_html: true
output:
  word_document:
    fig_caption: yes
    fig_width: 4
    keep_md: yes
  header-includes:
  - \usepackage{knputenc}
  - \usepackage{unicode-math}
  - \pagenumbering{gobble}
  - \documentclass{article}
  - \usepackage{amsmath}
  pdf_document:
    toc: FALSE
    toc_depth: 3
    fig_caption: yes
    fig_height: 4
    fig_width: 5
    highlight: tango
    keep_tex: yes
    latex_engine: xelatex
  html_document:
    df_print: paged
    toc: true
    toc_depth: 3
    number_sections: true
    toc_float:
      collapsed: false
      smooth_scroll: true
#runtime: shiny
---

```{r startup, eval=TRUE, echo=FALSE, results='hide',message=FALSE}
 
 #source("R/make.R")       # loads packages, data, setup, etc.
 knitr::opts_chunk$set(echo = T, fig.align="center",warning = FALSE, message = FALSE)
 #knitr::opts_knit$set(root.dir = '../')
 #options(knitr.table.format = "latex")
 thisYr <- format(Sys.time(), "%Y")
 today  <- format(Sys.time(), "%b %d, %Y")
 
suppressWarnings(source("R/make.R"))
load("ESR_EBS/endofCent.Rdata")
load("ESR_EBS/midCent.Rdata")
load("ESR_EBS/endofCent_N.Rdata")
load("ESR_EBS/midCent_N.Rdata")

 BT  <- endofCent%>%filter(lab2=="BT")
 SST <- endofCent%>%filter(lab2=="SST")
  
 BTN  <- endofCentN%>%filter(lab2=="BT")
 SSTN <- endofCentN%>%filter(lab2=="SST")
 
```

```{r echo=FALSE, eval=FALSE, include = FALSE,results='hide'}

    rm(list=ls())
    # load data, packages, and functions
    # ------------------------------------
    source("R/make.R")

    # use this to get vals for the table for current_yr+1 and current_yr +2
    thisYr <- as.numeric(format(Sys.time(), "%Y"))
    today  <- format(Sys.time(), "%b %d, %Y")
    
    # setwd("D:/GitHub_cloud/ACLIM2")
    # The three CMIP5 models used as the basis for this study are:
    # GFDL-ESM2M (Dunne et al., 2012), CESM (Kay et al., 2015), and MIROC (Watanabe et al., 2011).
    # Output from three different global Earth System Models 
    # (ESMs) were used in this study: 1) CESM version 2 with 
    # Community Atmospheric Model version 6 (CESM2-CAM6; Danabasoglu et al., 2020); 2) 
    # GFDL Earth System Model version 4.1 (GFDL-ESM4; Dunne et al., 2020); 3) 
    # MIROC-Earth System version 2 for Long-term simulations (MIROC-ES2L; Hajima et al., 2020). 
    
    # Create lookup table:
    # ------------------------
    
     
       cmip6_lkup <- data.frame(rbind(
       c( "MIROC-ES2L", " r1i1p1f2","miroc"),
       c("GFDL-ESM4", " r1i1p1f1","gfdl"),
       c("CESM2", " r2i1p1f1","cesm")),stringsAsFactors = FALSE)
       colnames(cmip6_lkup)<-c("model", "ensemble","GCM")
      

      cmip5_lkup <- data.frame(rbind(
       c( "MIROC5", " r1i1p1","miroc"),
       c("GFDL-ESM2M", " r1i1p1","gfdl"),
       c("CESM1-BGC", " r1i1p1","cesm")),stringsAsFactors = FALSE)
      colnames(cmip5_lkup)<-c("model", "ensemble","GCM")
      
      
    
    # Now get CMIP6 warming levels from Mathis repo
    # ------------------------
    # https://github.com/mathause/cmip_warming_levels
    
     GWL_cmip6 <- read.csv("Data/in/cmip_warming_levels-v0.2.0/mathause-cmip_warming_levels-fc9a5d7/warming_levels/cmip6_all_ens/csv/cmip6_warming_levels_all_ens_1850_1900.csv",skip=4)
     GWL_cmip5 <- read.csv("Data/in/cmip_warming_levels-v0.2.0/mathause-cmip_warming_levels-fc9a5d7/warming_levels/cmip5_all_ens/csv/cmip5_warming_levels_all_ens_1850_1900.csv",skip=4)
     GWL_cmip6%>%filter(model%in%c("CESM2","MIROC-ES2L","GFDL-ESM4"),exp%in%c(" ssp126"," ssp585"))
     GWL_cmip5%>%filter(model%in%c("CESM1-BGC","MIROC5","GFDL-ESM4"),exp%in%c(" rcp45"," rcp85"))
    
      
      cmip5<- GWL_cmip5%>%
        filter(model%in%cmip5_lkup$model,exp%in%c(" rcp45"," rcp85"))%>%
        left_join(cmip5_lkup,by=c("model"="model","ensemble"="ensemble"))%>%
        filter(!is.na(GCM))%>%
        group_by(model,exp,warming_level,GCM)%>%summarize(
          mnSTyr = mean(start_year),
          mnENDyr = mean(end_year))%>%mutate(exp = gsub(" ","",exp))
      
      cmip6 <- GWL_cmip6%>%filter(
        model%in%c("CESM2","MIROC-ES2L","GFDL-ESM4"),
        exp%in%c(" ssp126"," ssp585"))%>%
           left_join(cmip6_lkup,by=c("model"="model","ensemble"="ensemble"))%>%
         filter(!is.na(GCM))%>%
        group_by(model,exp,warming_level,GCM)%>%summarize(
          mnSTyr = mean(start_year),
          mnENDyr = mean(end_year))%>%mutate(exp = gsub(" ","",exp))
       
    # Load the rest of the ROMSNPZ data:
    # ------------------------
    cat("Rdata_path : ", Rdata_path,"\n") 
    
    hind_yearsIN    <- 1979:thisYr
    proj_yearsIN    <- (thisYr+1):2100

    yrs             <- hind_yearsIN
    plot_vars      <- c("fracbelow2",
                        "fracbelow0",
                        "temp_bottom5m",
                        "temp_surface5m",
                        "oxygen_bottom5m",
                        "pH_bottom5m",
                        "pH_surface5m",
                        "aice","EupS_integrated","NCaS_integrated")

    # Load weekly data:
    # ----------------------------------------------
    load("Data/out/K20P19_CMIP5/allEBS_means/ACLIM_weekly_fut_mn.Rdata")
      ACLIM_weekly_futCMIP5 <- ACLIM_weekly_fut%>%left_join(weekly_var_def,join_by(var== name)) 
    load("Data/out/K20P19_CMIP6/allEBS_means/ACLIM_weekly_fut_mn.Rdata")
    load("Data/out/K20P19_CMIP6/allEBS_means/ACLIM_weekly_hind_mn.Rdata")
    load("Data/out/K20P19_CMIP6/allEBS_means/ACLIM_weekly_hist_mn.Rdata")
    ACLIM_weekly_hist         <- ACLIM_weekly_hist%>%left_join(weekly_var_def,join_by(var== name)) 
    ACLIM_weekly_fut          <- ACLIM_weekly_fut%>%left_join(weekly_var_def,join_by(var== name)) 
    ACLIM_weekly_futCMIP6     <- ACLIM_weekly_fut
    ACLIM_weekly_futCMIP5$GCM <- gsub("MIROC","miroc",ACLIM_weekly_futCMIP5$GCM)
    ACLIM_weekly_futCMIP5$GCM <- gsub("CESM","cesm",  ACLIM_weekly_futCMIP5$GCM)
    ACLIM_weekly_futCMIP5$GCM <- gsub("GFDL","gfdl",  ACLIM_weekly_futCMIP5$GCM)
    ACLIM_weekly_fut          <- rbind(ACLIM_weekly_futCMIP6,ACLIM_weekly_futCMIP5)
    
    # Load surveyrep data:
    # ----------------------------------------------
    load("Data/out/K20P19_CMIP5/allEBS_means/ACLIM_surveyrep_fut_mn.Rdata")
    ACLIM_surveyrep_futCMIP5     <- ACLIM_surveyrep_fut%>%left_join(srvy_var_def,join_by(var== name)) 
    load("Data/out/K20P19_CMIP6/allEBS_means/ACLIM_surveyrep_fut_mn.Rdata")
    load("Data/out/K20P19_CMIP6/allEBS_means/ACLIM_surveyrep_hind_mn.Rdata")
    load("Data/out/K20P19_CMIP6/allEBS_means/ACLIM_surveyrep_hist_mn.Rdata")
    ACLIM_surveyrep_hist         <- ACLIM_surveyrep_hist%>%left_join(srvy_var_def,join_by(var== name)) 
    ACLIM_surveyrep_hind         <- ACLIM_surveyrep_hind%>%left_join(srvy_var_def,join_by(var== name)) 
    ACLIM_surveyrep_fut          <- ACLIM_surveyrep_fut%>%left_join(srvy_var_def,join_by(var== name)) 
    ACLIM_surveyrep_futCMIP6     <- ACLIM_surveyrep_fut
    ACLIM_surveyrep_futCMIP5$GCM <- gsub("MIROC","miroc",ACLIM_surveyrep_futCMIP5$GCM)
    ACLIM_surveyrep_futCMIP5$GCM <- gsub("CESM","cesm",  ACLIM_surveyrep_futCMIP5$GCM)
    ACLIM_surveyrep_futCMIP5$GCM <- gsub("GFDL","gfdl",  ACLIM_surveyrep_futCMIP5$GCM)
    ACLIM_surveyrep_fut          <- rbind(ACLIM_surveyrep_futCMIP6,ACLIM_surveyrep_futCMIP5)
     
    # Load annual data:
    # ----------------------------------------------    
    load("Data/out/K20P19_CMIP5/allEBS_means/ACLIM_annual_fut_mn.Rdata")
    ACLIM_annual_futCMIP5     <- ACLIM_annual_fut%>%left_join(weekly_var_def,join_by(var== name)) 
    load("Data/out/K20P19_CMIP6/allEBS_means/ACLIM_annual_fut_mn.Rdata")
    load("Data/out/K20P19_CMIP6/allEBS_means/ACLIM_annual_hind_mn.Rdata")
    load("Data/out/K20P19_CMIP6/allEBS_means/ACLIM_annual_hist_mn.Rdata")
    ACLIM_annual_hist         <- ACLIM_annual_hist%>%left_join(weekly_var_def,join_by(var== name)) 
    ACLIM_annual_hind         <- ACLIM_annual_hind%>%left_join(weekly_var_def,join_by(var== name)) 
    ACLIM_annual_fut          <- ACLIM_annual_fut%>%left_join(weekly_var_def,join_by(var== name)) 
    ACLIM_annual_futCMIP6     <- ACLIM_annual_fut
    ACLIM_annual_futCMIP5$GCM <- gsub("MIROC","miroc",ACLIM_annual_futCMIP5$GCM)
    ACLIM_annual_futCMIP5$GCM <- gsub("CESM","cesm",  ACLIM_annual_futCMIP5$GCM)
    ACLIM_annual_futCMIP5$GCM <- gsub("GFDL","gfdl",  ACLIM_annual_futCMIP5$GCM)
    ACLIM_annual_fut          <- rbind(ACLIM_annual_futCMIP6,ACLIM_annual_futCMIP5)
      
    # Load seasonal data:
    # ----------------------------------------------    
    load("Data/out/K20P19_CMIP5/allEBS_means/ACLIM_seasonal_fut_mn.Rdata")
    ACLIM_seasonal_futCMIP5     <- ACLIM_seasonal_fut%>%left_join(weekly_var_def,join_by(var== name)) 
    load("Data/out/K20P19_CMIP6/allEBS_means/ACLIM_seasonal_fut_mn.Rdata")
    load("Data/out/K20P19_CMIP6/allEBS_means/ACLIM_seasonal_hind_mn.Rdata")
    load("Data/out/K20P19_CMIP6/allEBS_means/ACLIM_seasonal_hist_mn.Rdata")
    ACLIM_seasonal_hist         <- ACLIM_seasonal_hist%>%left_join(weekly_var_def,join_by(var== name)) 
    ACLIM_seasonal_hind         <- ACLIM_seasonal_hind%>%left_join(weekly_var_def,join_by(var== name)) 
    ACLIM_seasonal_fut          <- ACLIM_seasonal_fut%>%left_join(weekly_var_def,join_by(var== name)) 
    ACLIM_seasonal_futCMIP6     <- ACLIM_seasonal_fut
    ACLIM_seasonal_futCMIP5$GCM <- gsub("MIROC","miroc",ACLIM_seasonal_futCMIP5$GCM)
    ACLIM_seasonal_futCMIP5$GCM <- gsub("CESM","cesm",  ACLIM_seasonal_futCMIP5$GCM)
    ACLIM_seasonal_futCMIP5$GCM <- gsub("GFDL","gfdl",  ACLIM_seasonal_futCMIP5$GCM)
    ACLIM_seasonal_fut          <- rbind(ACLIM_seasonal_futCMIP6,ACLIM_seasonal_futCMIP5)
    
    # Load monthly data:
    # ----------------------------------------------
    load("Data/out/K20P19_CMIP5/allEBS_means/ACLIM_monthly_fut_mn.Rdata")
    ACLIM_monthly_futCMIP5 <- ACLIM_monthly_fut%>%left_join(weekly_var_def,join_by(var== name)) 
    load("Data/out/K20P19_CMIP6/allEBS_means/ACLIM_monthly_fut_mn.Rdata")
    load("Data/out/K20P19_CMIP6/allEBS_means/ACLIM_monthly_hind_mn.Rdata")
    load("Data/out/K20P19_CMIP6/allEBS_means/ACLIM_monthly_hist_mn.Rdata")
    ACLIM_monthly_hist         <- ACLIM_monthly_hist%>%left_join(weekly_var_def,join_by(var== name)) 
      ACLIM_monthly_hind         <- ACLIM_monthly_hind%>%left_join(weekly_var_def,join_by(var== name))
    ACLIM_monthly_fut          <- ACLIM_monthly_fut%>%left_join(weekly_var_def,join_by(var== name)) 
    ACLIM_monthly_futCMIP6     <- ACLIM_monthly_fut
    ACLIM_monthly_futCMIP5$GCM <- gsub("MIROC","miroc",ACLIM_monthly_futCMIP5$GCM)
    ACLIM_monthly_futCMIP5$GCM <- gsub("CESM","cesm",  ACLIM_monthly_futCMIP5$GCM)
    ACLIM_monthly_futCMIP5$GCM <- gsub("GFDL","gfdl",  ACLIM_monthly_futCMIP5$GCM)
    ACLIM_monthly_fut          <- rbind(ACLIM_monthly_futCMIP6,ACLIM_monthly_futCMIP5)
   
    #ACLIM_weekly_fut<-rbind(ACLIM_weekly_fut,ACLIM_weekly_futCMIP5,ACLIM_weekly_futCMIP5H16)
   
    allyears   <- min(ACLIM_weekly_hind$year,na.rm=T):max(ACLIM_weekly_fut$year,na.rm=T)
    hindyrs    <- min(ACLIM_weekly_hind$year,na.rm=T):max(ACLIM_weekly_hind$year,na.rm=T)
    futyrs     <- min(ACLIM_weekly_fut$year,na.rm=T):max(ACLIM_weekly_fut$year,na.rm=T)
    scens      <- c("ssp126"  ,"rcp45"  ,"ssp585","rcp85")
    GCMs       <- c("miroc", "gfdl" , "cesm" )
    
    
  # Create climatology
  # -------------------------------------------
  # get just the core columns:
  grplist <- c("var","sim",
               "basin","type",
               "sim_type",
               "gcmcmip","CMIP",
               "scen",
               "GCM","mod")
  
  plot_vars      <- c("fracbelow2",
                      "fracbelow0",
                      "temp_bottom5m",
                      "temp_surface5m",
                      "oxygen_bottom5m",
                      "pH_bottom5m",
                      "pH_surface5m",
                      "aice","EupS_integrated","NCaS_integrated")

  hind <- ACLIM_weekly_hind%>%filter(var%in%plot_vars)%>%
    select(all_of(c(grplist,"year","season","mn_val","mnVal_hind","mnDate","jday","wk","mo")))
  
  # get weekly climatology
  wkly_clim <- hind%>%ungroup()%>%filter(year<2019)%>%
    group_by(!!!syms(c(grplist,"wk")))%>%
    summarise(
      jday_clim = mean(jday, na.rm=T),
      mn_clim = mean(mn_val,na.rm=T),
      sd_clim = sd(mn_val,na.rm=T))%>%
    mutate(mndate_clim = as.Date("1990-01-01")+jday_clim)
  ggplot(wkly_clim%>%filter(var=="EupS_integrated"))+
    geom_line(aes(x=mndate_clim,y= mn_clim,color=basin))+
    geom_line(aes(x=mndate_clim,y= mn_clim+sd_clim,color=basin),linetype="dashed")+
    geom_line(aes(x=mndate_clim,y= mn_clim-sd_clim,color=basin),linetype="dashed")
  
  
  # get monthly climatology
  mo_clim <- hind%>%ungroup()%>%filter(year<2019)%>%
    group_by(!!!syms(c(grplist,"mo")))%>%
    summarise(
      jday_clim = mean(jday, na.rm=T),
      mn_clim = mean(mn_val,na.rm=T),
      sd_clim = sd(mn_val,na.rm=T))%>%
    mutate(mndate_clim = as.Date("1990-01-01")+jday_clim)
  ggplot(mo_clim%>%filter(var=="EupS_integrated"))+
    geom_line(aes(x=mndate_clim,y= mn_clim,color=basin))+
    geom_line(aes(x=mndate_clim,y= mn_clim+sd_clim,color=basin),linetype="dashed")+
    geom_line(aes(x=mndate_clim,y= mn_clim-sd_clim,color=basin),linetype="dashed")
  
  
  # get seasonal climatology
  season_clim <- hind%>%ungroup()%>%filter(year<2019)%>%
    group_by(!!!syms(c(grplist,"season")))%>%
    summarise(
      jday_clim = mean(jday, na.rm=T),
      mn_clim = mean(mn_val,na.rm=T),
      sd_clim = sd(mn_val,na.rm=T))%>%
    mutate(mndate_clim = as.Date("1990-01-01")+jday_clim)
  ggplot(season_clim%>%filter(var=="EupS_integrated"))+
    geom_line(aes(x=mndate_clim,y= mn_clim,color=basin))+
    geom_line(aes(x=mndate_clim,y= mn_clim+sd_clim,color=basin),linetype="dashed")+
    geom_line(aes(x=mndate_clim,y= mn_clim-sd_clim,color=basin),linetype="dashed")
  
  # get annual climatatology
  annual_clim <- hind%>%ungroup()%>%filter(year<2019)%>%
    group_by(!!!syms(c(grplist)))%>%
    summarise(
      jday_clim = mean(jday, na.rm=T),
      mn_clim = mean(mn_val,na.rm=T),
      sd_clim = sd(mn_val,na.rm=T))%>%
    mutate(mndate_clim = as.Date("1990-01-01")+jday_clim)

    
    # Create timeblocks
    # -------------------------------------------
    blocks      <- round(unique(round(allyears/1000,2)*1000))
    if(any(blocks>max(futyrs))) 
      blocks    <- blocks[-which(blocks>max(futyrs))]
    blocks_fut  <- unique(round(futyrs/1000,2)*1000)
    blocks_hind <- unique(10+round(hindyrs/1000,2)*1000)
    subyear     <- blocks
    
    grplist <- c("var","sim",
                 "basin","type",
                 "sim_type","jday",
                 "gcmcmip","CMIP",
                 "GCM","mod")
    
    length.na <- function(x,na.rm=F){
      if(na.rm){
        x <- as.numeric(na.omit(x))
      }
      length(x)
    }
    
    se <- function(x, na.rm=T){
      if(na.rm){
        x <- as.numeric(na.omit(x))
      }
      if(length(x)>0){
        sd(x)/sqrt(length(x))
      }else{
        0
      }
      
    }
    
    # Re-org the data for plots
    # ----------------------------------------------
    dat <- 
      ACLIM_weekly_fut%>%
      filter(var%in%plot_vars)%>%
      mutate(yrFactor   = factor(as.character(year),levels= allyears),
             year_block = 10+round(year/1000,2)*1000)%>%
      mutate(year_block = factor(as.character(year_block),levels=blocks),
              scen      = factor(RCP,levels= scens))
    
    
    dat_hind   <- ACLIM_weekly_hind%>%
      filter(var%in%plot_vars)%>%
      mutate(yrFactor = factor(as.character(year),levels= allyears),
             year_block= 10+round(year/1000,2)*1000)%>%
      mutate(year_block = factor(as.character(year_block),levels=blocks))
    
    dat_hist   <- ACLIM_weekly_hist%>%
      filter(var%in%plot_vars)%>%
      mutate(yrFactor = factor(as.character(year),levels= allyears),
             year_block= 10+round(year/1000,2)*1000)%>%
      mutate(year_block = factor(as.character(year_block),levels=blocks),
             scen      = factor(RCP,levels= scens))
    
    mndat<-
      dat_hind%>%group_by(!!!syms(grplist))%>%
      summarize(mnVal_hind2 = mean(mn_val,na.rm = T),
                mn_hind    = mean(mnVal_hind,na.rm = T),
                sdVal_hind = sd(mn_val,na.rm = T),
                seVal_hind = se(mn_val,na.rm = T),
                n = length.na(mn_val,na.rm = T))
     histdat<-
      dat_hist%>%group_by(!!!syms(grplist))%>%
      summarize(mnVal_hist2 = mean(mn_val,na.rm = T),
                mn_hist     = mean(mnVal_hist,na.rm = T),
                sd_val_hist = sd(mn_val,na.rm = T),
                se_val_hist = se(mn_val,na.rm = T),
                n_hist      = length.na(mn_val,na.rm = T))
     
     
    
     dat_hind <- dat_hind%>%left_join(wkly_clim)
      dat_hind <- dat_hind%>%left_join(mndat)%>%ungroup()%>%data.frame()
     
      #get mean vals by timeblock
    dat_b <- dat%>%
      group_by(basin ,GCM, sim,scen, 
               var,jday,
               longname,units,
               mnVal_hind,year_block )%>%
      summarize(
        mn_val            = mean(mn_val,na.rm=T),
        val_biascorrected = mean(val_biascorrected,na.rm=T))%>%
      rename(yrFactor     = year_block)
    
    dat_hind_b <- dat_hind%>%
      group_by(basin ,GCM, sim,scen, 
               var,jday,
               sdVal_hind,
               mnVal_hind,
               year_block )%>%
      summarize(
        mn_val = mean(mn_val,na.rm=T))%>%
      rename(yrFactor = year_block)
    
    
    dat_s      <- ACLIM_seasonal_fut%>%
      mutate(scen = RCP,
      gcmscen     = paste0(GCM,"_",scen),
      season_var  = paste0(season,"_",var))
    
    ACLIM_seasonal_hind <- ACLIM_seasonal_hind%>%
      mutate(
      gcmscen     = paste0(GCM,"_",scen),
      season_var  = paste0(season,"_",var))%>%
      group_by(basin ,GCM, sim,scen, 
               var,season,
               longname,units,
               mnVal_hind)%>%
      summarize(
        mnVal_hind2 = mean(mn_val,na.rm=T),
        sdVal_hind = sd(mn_val,na.rm=T))
    
    ACLIM_monthly_hind <- ACLIM_monthly_hind%>%
      mutate(
      gcmscen     = paste0(GCM,"_",scen),
      season_var  = paste0(season,"_",var))%>%
      group_by(basin ,GCM, sim,scen, 
               var,mo,
               longname,units,
               mnVal_hind)%>%
      summarize(
        mnVal_hind2 = mean(mn_val,na.rm=T),
        sdVal_hind  = sd(mn_val,na.rm=T))
    
    dat_m      <- ACLIM_monthly_fut%>%
      mutate(scen=RCP,
      gcmscen = paste0(GCM,"_",scen),
      mo_var  = paste0(mo,"_",var))
   
    i<- 0
    for(s in 1:length(scens)){
      for(g in 1:length(GCMs)){
        i <- i +1
        if(i == 1){
          dat_mn_hind    <- dat_hind%>%
            mutate(scen = factor(scens[s],levels= scens), GCM = GCMs[g])
          dat_mn_hind_b  <- dat_hind_b%>%
            mutate(scen = factor(scens[s],levels= scens), GCM = GCMs[g])
          dat_mn_s_hind  <- ACLIM_seasonal_hind%>%
            mutate(scen = factor(scens[s],levels= scens), GCM = GCMs[g])
          dat_mn_m_hind  <- ACLIM_monthly_hind%>%
            mutate(scen = factor(scens[s],levels= scens), GCM = GCMs[g])
        
        }else{
          dat_mn_hind    <- rbind(dat_mn_hind,dat_hind%>%
                                    mutate(scen = factor(scens[s],levels= scens), GCM = GCMs[g]))
          dat_mn_hind_b  <- rbind(dat_mn_hind_b,dat_hind_b%>%
                                    mutate(scen = factor(scens[s],levels= scens), GCM = GCMs[g]))
          dat_mn_s_hind  <- rbind(dat_mn_s_hind,ACLIM_seasonal_hind%>%
                                    mutate(scen = factor(scens[s],levels= scens), GCM = GCMs[g]))
          dat_mn_m_hind  <- rbind(dat_mn_m_hind,ACLIM_monthly_hind%>%
                                    mutate(scen = factor(scens[s],levels= scens), GCM = GCMs[g]))
        }
      }
    }
    
    dat_mn_s_hind <- dat_mn_s_hind%>%mutate(
      gcmscen    = paste0(GCM,"_",scen),
      season_var = paste0(season,"_",var))
    
    dat_mn_m_hind <- dat_mn_m_hind%>%mutate(
      gcmscen = paste0(GCM,"_",scen),
      mo_var  = paste0(mo,"_",var))
    
    dat_mn_hind_sub <- dat_mn_hind%>%filter(year == 2010)%>%
      mutate(yrFactor = factor("[1980-2010] mean",levels=c("[1980-2010] mean")))
    
    dat_sub<- dat%>%filter(year%in%subyear, year>2010)%>%
      mutate(yrFactor = factor(as.character(year),levels =subyear))
    
      scensIN      <- c("ssp126","ssp585")
      GCMsIN       <- GCMs
      thinline     <- 0.7
      alpha_annual <- 0.6
      alpha_blocks <- 1

     pvar2         <- c("temp_bottom5m","fracbelow2","pH_bottom5m")
     pseason       <- c("Summer","Summer","Winter")
    
     pvar          <- c("temp_bottom5m","pH_bottom5m")
     pseason       <- c("Summer","Winter")
     
     pvar2         <- c("temp_surface5m","temp_bottom5m")
     pseason2      <- c("Summer","Summer")
     
     pvar3         <- c("temp_surface5m","temp_bottom5m","NCaS_integrated",
                        "NCaS_integrated","EupS_integrated","aice","fracbelow2")
     pseason3      <- c("Summer","Summer","Spring","Fall","Fall","Spring","Summer")
     
     season_varIN  <- paste0(pseason,"_",pvar)
     season_varIN  <- paste0(pseason2,"_",pvar2)
      names(season_varIN) <- c(paste(pseason2,c("SST", "BT")))
     season_varIN  <- paste0(pseason3,"_",pvar3)
      names(season_varIN) <- c(paste(pseason2,c("SST", "BT","Neocalanus",
                        "Neocalanus","Euphausiid","Sea ice","Coldpool")))
    
     
     
     ylablkup     <- melt(season_varIN)%>%rename(season_var = value)
     ylablkup$lab <- rownames(ylablkup)
     
     ylablkup <- data.frame(var = 
                              c("temp_surface5m","temp_surface5m",
                                "temp_bottom5m","temp_bottom5m"),
                            season = c("Winter","Summer","Winter","Summer"),
                            varlab = c("SST","SST","BT","BT"))
     ylablkup<- ylablkup%>%mutate(season_var = paste0(season,"_",var), lab = paste(season, varlab))
     
    # Now make climatology plots
    # ------------------------------------
     pp_BT <-   plot_weekly(   datIN = dat, 
                               dat_hindIN = dat_mn_hind, 
                               datIN_sub  = dat_sub,
                               datIN_mn_hind_sub = dat_mn_hind_sub,
                               useBiascorrected=FALSE,
                               varIN = "temp_bottom5m",  basinIN = "SEBS",
                               titleIN =  "SEBS Bottom temperature")
     pp_pH <-   plot_weekly( varIN = "pH_bottom5m",    basinIN = "SEBS",
                             titleIN =  "SEBS Bottom pH")
     pp_O2 <-   plot_weekly( varIN = "oxygen_bottom5m",basinIN = "SEBS",
                             titleIN =  "SEBS Bottom O2")
     pp_CP <-   plot_weekly( varIN = "fracbelow2",     basinIN = "SEBS",
                             titleIN =  "SEBS Cold pool proxy ( fraction < 2 oC)")
    
     pp_BTN <-   plot_weekly( varIN = "temp_bottom5m",  basinIN = "NEBS",
                             titleIN =  "NEBS Bottom temperature")
     pp_pHN <-   plot_weekly( varIN = "pH_bottom5m",    basinIN = "NEBS",
                             titleIN =  "NEBS Bottom pH")
     pp_O2N <-   plot_weekly( varIN = "oxygen_bottom5m",basinIN = "NEBS",
                             titleIN =  "NEBS Bottom O2")
     pp_CPN <-   plot_weekly( varIN = "fracbelow2",     basinIN = "NEBS",
                             titleIN =  "NEBS Cold pool proxy ( fraction < 2 oC)")
     
     
     # Now make annual TS
     # ------------------------------------
     season_var1 <- c('Summer SST (oC)' ="Summer_temp_surface5m",
                      'Summer BT (oC)' = "Summer_temp_bottom5m" )
     season_var2 <- c('Summer BT (oC)' ="Summer_temp_bottom5m",
                      'Winter pH (bottom)' = "Winter_pH_bottom5m" )
     
     season_var3 <- c('Summer BT (oC)'    ="Summer_temp_bottom5m",
                      'Summer coldpool'    ="Summer_fracbelow2",
                      'Winter pH (bottom)' = "Winter_pH_bottom5m",
                      'Winter sea ice' = "Winter_aice"
                      )
     
   
      test  <- plot_annual(  datIN  = dat_s,
                             datIN_hind = dat_mn_s_hind,
                             basinIN = "NEBS",
                             season_varIN = c('Summer SST' ="Summer_temp_surface5m"), 
                             futline =.5,alpha_fut =.5,alpha_loess = .4)
     pp_NEBS2  <- plot_annual(  datIN  = dat_s,
                                datIN_hind = dat_mn_s_hind,
                              basinIN = "NEBS",season_varIN = season_var2, 
                              futline =.5,alpha_fut =.5,alpha_loess = .4)
     pp_SEBS2  <- plot_annual(basinIN = "SEBS",season_varIN = season_var2)
     pp_NEBS   <- plot_annual(basinIN = "NEBS",season_varIN = season_var1)
     pp_SEBS   <- plot_annual(basinIN = "SEBS",season_varIN = season_var1)
     
     pp_NEBS3   <- plot_annual(basinIN = "NEBS",season_varIN = season_var3)
     pp_SEBS3   <- plot_annual(basinIN = "SEBS",season_varIN = season_var3)
   
     # Now make plot of warm~GWL
     # ------------------------------------
      gwls <- make_GWLplot(ylablkupIN = ylablkup,
                           cmipIN = cmip6, 
                           datIN  = dat_s,
                           datIN_hind = dat_mn_s_hind)
       	 
     varlab2 <- data.frame(var= c("temp_surface5m","temp_bottom5m"),
                          lab2 = c("SST","BT"))
     
     endofCentN <- pp_NEBS$fut%>%ungroup()%>%filter(year>2079)%>%
       left_join(varlab2)%>%
       group_by(scen,basin, var,lab2)%>%
        summarize(mn = mean(val_biascorrected),
                  min =min(val_biascorrected),
     	            max=max(val_biascorrected),
                  sd = sd(val_biascorrected), 
                  hind = mean(mnVal_hind),
                  hindsd = mean(sdVal_hind),
                  hist = mean(mnVal_hist),
                  mndelta = mean(mn_val-mnVal_hist), 
                  sddelta = sd(mn_val-mnVal_hist),
                  mindelta = min(mn_val-mnVal_hist),
                  maxdelta = max(mn_val-mnVal_hist))
     
     midCentN   <- pp_NEBS$fut%>%ungroup()%>%filter(year>2039&year<2061)%>%
        left_join(varlab2)%>%
       group_by(scen,basin, var,lab2)%>%
        summarize(mn  = mean(val_biascorrected), 
                  min = min(val_biascorrected),
     	            max = max(val_biascorrected),
                  sd  = sd(val_biascorrected), 
                  hind = mean(mnVal_hind),
                  hindsd = mean(sdVal_hind),
                  hist = mean(mnVal_hist),
                  mndelta = mean(mn_val-mnVal_hist), 
                  sddelta = sd(mn_val-mnVal_hist),
                  mindelta = min(mn_val-mnVal_hist),
                  maxdelta = max(mn_val-mnVal_hist))
     
     endofCent  <- pp_SEBS$fut%>%ungroup()%>%filter(year>2079)%>%
        left_join(varlab2)%>%
       group_by(scen,basin, var,lab2)%>%
     	  summarize(mn = mean(val_biascorrected), 
     	            min =min(val_biascorrected),
     	            max=max(val_biascorrected),
     	            sd = sd(val_biascorrected), 
     	            hind = mean(mnVal_hind),
                  hindsd = mean(sdVal_hind),
     	            hist = mean(mnVal_hist),
                  mndelta = mean(mn_val-mnVal_hist), 
                  sddelta = sd(mn_val-mnVal_hist),
                  mindelta = min(mn_val-mnVal_hist),
                  maxdelta = max(mn_val-mnVal_hist))
     
     midCent    <- pp_SEBS$fut%>%ungroup()%>%filter(year>2039&year<2061)%>%
     	   left_join(varlab2)%>%
       group_by(scen,basin, var,lab2)%>%
       summarize(mn = mean(val_biascorrected), 
                 min =min(val_biascorrected),
     	            max=max(val_biascorrected),
                 sd = sd(val_biascorrected), 
                  hind = mean(mnVal_hind),
                  hindsd = mean(sdVal_hind),
                 hist = mean(mnVal_hist),
                  mndelta = mean(mn_val-mnVal_hist),
                 sddelta = sd(mn_val-mnVal_hist),
                  mindelta = min(mn_val-mnVal_hist),
                 maxdelta = max(mn_val-mnVal_hist))
     
     flout <- file.path("ESR_EBS",as.Date(now()))
     if(!dir.exists(flout))
       dir.create(flout)
     
     if(update.figs){
     
      save(endofCent, file = file.path(flout,"endofCent.Rdata"))
      save(midCent,   file = file.path(flout,"midCent.Rdata"))
      save(endofCentN,file = file.path(flout,"endofCent_N.Rdata"))
      save(midCentN,  file = file.path(flout,"midCent_N.Rdata"))
      
      # Make figs directory:
      if(!dir.exists(file.path(flout,"Figs"))) dir.create(file.path(flout,"Figs"))
      
      # save plots:
      # -------------------------------
     	sclr <-1
      save_png(pp_NEBS2$pp,  file = file.path(flout,"annualTS_NEBS.png"),w = 6*sclr, h = 6*sclr)
      save_png(pp_NEBS$pp, file = file.path(flout,"ACLIM_Holsman_Fig2.png"),w = 6*sclr, h = 6*sclr)
      save_png(pp_SEBS2$pp,  file= file.path(flout,"annualTS.png"),w = 6*sclr, h = 6*sclr)
      save_png(pp_SEBS$pp, file = = file.path(flout,"ACLIM_Holsman_Fig1.png"),w = 6*sclr, h = 6*sclr)
     
      sclr <-1
      save_png( gwls$pp, file = file.path(flout,"ACLIM_Holsman_Fig3.png"),w = 6*sclr, h = 6*sclr)
     
      sclr <-1.3
      save_png( pp_BT, file =  file.path(flout,"ACLIM_Holsman_Fig4.png"),w = 6*sclr, h = 6*sclr)
      save_png( pp_pH, file = file.path(flout,"nonBC_weeklyProj_pH.png"),w = 6*sclr, h = 6*sclr)
      save_png( pp_O2, file = file.path(flout,"nonBC_weeklyProj_O2.png"),w = 6*sclr, h = 6*sclr)
      save_png( pp_CP, file = "ESR_EBS/Figs/nonBC_weeklyProj_CP.png",w = 6*sclr, h = 6*sclr)
      
      save_png( pp_BTN, file  = file.path(flout,"ACLIM_Holsman_Fig5.png"),w = 6*sclr, h = 6*sclr)
      save_png( pp_pHN, file  = file.path(flout,"nonBC_weeklyProj_pH_N.png"),w = 6*sclr, h = 6*sclr)
      save_png( pp_O2N, file  = file.path(flout,"nonBC_weeklyProj_O2_N.png"),w = 6*sclr, h = 6*sclr)
      save_png( pp_CPN, file  = file.path(flout,"nonBC_weeklyProj_CP_N.png"),w = 6*sclr, h = 6*sclr)
      }
      
  #     
  # 	# now set up the docs folder
  # 	EcoCons_path <- file.path("docs/EcoCons")
  # 	data_oath  <- file.path("docs/MSM_Assessment")
  # 	
  # 	if(!dir.exists(file.path(proj_dir,proj,"docs")))
  # 	   dir.create(file.path(proj_dir,proj,"docs"))
  # 	if(!dir.exists(EcoCons_path))
  # 	   dir.create(EcoCons_path)
  # 	if(!dir.exists(assmnt_path))
  # 	   dir.create(assmnt_path)
  # 	
  # 	main       <- (file.path(proj_dir,proj,"runs",fldr_use))
  # 	data_fl    <- file.path(proj_dir,proj,"data/in/dat_input_files")
  # 	fl_0       <- file.path(main,paste0(filenm,"_0"))
  # 	fl_2       <- file.path(main,paste0(filenm,"_2"))
  #   load(file.path(fl_0,"R_figures/plotdata.Rdata"))
  #   load(file.path(fl_2,"results/Mort.Rdata"))
  #   
  #   dir("Data/out/K20P19_CMIP5/allEBS_means/ACLIM_weekly_fut_mn.Rdata")
  # 
  #   Ecodoc     <- EcoCons_path
  #   doc        <- assmnt_path
  #   
  #    if(!file.exists(file.path(Ecodoc,"Figs")))
  #       dir.create(file.path(Ecodoc,"Figs"))
  #    if(!file.exists(file.path(doc,"Figs")))
  #       dir.create(file.path(doc,"Figs"))
  #   setwd(Ecodoc)
    
```


## High resolution climate change projections for the Eastern Bering Sea

Kirstin K. Holsman$^1$,Albert Hermann$^2$, Wei Cheng$^2$, Kelly Kearney$^2$, Darren Pilcher$^2$, Kerim Aydin$^1$, Ivonne Ortiz$^2$

kirstin.holsman@noaa.gov

$^1$Cooperative Institute for Climate, Ocean and Ecosystem Studies, University of Washington, Seattle, WA. 98195

$^2$Alaska Fisheries Science Center, NOAA, 7600 Sand Point Way N.E., Bld. 4, Seattle, Washington 98115

**Last Updated: October `r thisYr`** 

## Summary statement:

A high resolution oceanographic model for the Bering Sea projects widespread warming across the region under low carbon mitigation (high greenhouse gas emission) scenarios ('ssp585' for Shared Socioeconomic Pathway 5-8.5) with modeled bottom temperatures exceeding historical ranges by mid-century onward. High carbon mitigation scenarios ('ssp126' for Shared Socioeconomic Pathway 1-2.6) are also associated with slight warming over the next century but modeled end-of-century temperatures do not exceed historical ranges (`r paste0(ref_years[1],"-", rev(ref_years)[1])`) as estimated by the hindcast. 

Global Warming Levels (GWL) of +3 and +4 $^o$C over pre-industrial temperatures (1850-1900) are associated with significant warming of surface and bottom water temperatures in both the Northern Bering and Southern Bering Seas (NEBS and SEBS, respectively). In contrast, under GWL of +1.5 $^o$C, SST and BT in both regions remain within the near-present range of climate variability (2010-2021). Global Warming Levels of +1.5 and +2 $^o$C represent the target and limit respectively of the UNFCCC [Paris Agreement and Nationally determined contributions (NDCs)](https://unfccc.int/ndc-synthesis-report-2022#Projected-GHG-Emission-levels); UNFCCC 2022) and are the warming levels beyond which climate change impacts and risks increase, --and feasibility and effectiveness of adaptation actions decrease--, rapidly with each GWL (IPCC 2022).

Warming of modeled bottom temperatures across seasons and models is projected in the SEBS under both mitigation scenarios, and in the NEBS under low carbon mitigation scenarios (ssp585) and is generally larger across regions and seasons under low carbon mitigation scenarios. While NEBS bottom water warming was also projected for summer months under high carbon mitigation scenarios (ssp126), only 1 of the 3 models projected substantial warming during winter months, indicating potential for sea ice and cold bottom water temperatures to be preserved to some extent through the end of century under high carbon mitigation scenarios (ssp126).

<!-- Under both low and high emission scenarios the model projects continued declines in ocean pH in bottom waters of the both the North and Southern Bering Sea, especially during winter months but declines are markedly larger and reach thresholds associated with decreased survival of crab larvae (7.8 Long et al. 20XX) by mid-century, and critical thresholds of high larval mortality by end of century (7.5 Long et al. 20xx).  -->

## Status and trends:  
 
Summer bottom temperatures in both the Southern Bering Sea (SEBS) and the Northern Bering Sea (NEBS) are projected to increase overtime, with higher rates of warming associated with low carbon mitigation (higher greenhouse gas emission) scenarios (ssp585) relative to high carbon mitigation scenarios (ssp126; Figs. 1, 2). Three Earth Systems Models (ESMs) are presented to reflect the spread in projections across ensemble members (see Description of Index below). There is general agreement in all three models (ESMs) with respect to trends in warming associated with alternative climate scenarios. For the SEBS, estimates of end of century warming of bottom temperatures [(2080-2100)-(`r paste0(ref_years[1],"-", rev(ref_years)[1])`)] range from `r a <- BT[BT$scen=="ssp126",]; paste0("+",round((a$mn-a$sd)-a$hind,2)," to +",round((a$mn+a$sd)-a$hind,2))` $^o$C for high carbon mitigation (ssp126) scenarios and `r a <- BT[BT$scen=="ssp585",]; paste0("+",round((a$mn-a$sd)-a$hind,2)," to +",round((a$mn+a$sd)-a$hind,2))` $^o$C for low carbon mitigation (ssp585) scenarios (ranges represent +/- 1 standard deviation; Fig. 1). For the NEBS, estimates of end of century warming of bottom temperatures [(2080-2100)-(`r paste0(ref_years[1],"-", rev(ref_years)[1])`)] range from `r a <- BTN[BTN$scen=="ssp126",]; paste0("+",round((a$mn-a$sd)-a$hind,2)," to +",round((a$mn+a$sd)-a$hind,2))` $^o$C for high carbon mitigation (ssp126) scenarios and `r a <- BTN[BTN$scen=="ssp585",]; paste0("+",round((a$mn-a$sd)-a$hind,2)," to +",round((a$mn+a$sd)-a$hind,2))` $^o$C for low carbon mitigation (ssp585) scenarios(Fig. 2). In high emission scenarios, bottom temperatures for the SEBS and NEBS by the mid-century (2050-2060) are projected to consistently (i.e., all three ensemble members) exceed the upper range of historical modeled temperatures from the Bering Sea shelf summer trawl survey and hindcast simulations. Mean historical (`r paste0(ref_years[1],"-", rev(ref_years)[1])`) bottom temperatures from the hindcast are `r round(BT$hind[1],1)` (SD = `r round(BT$hindsd[1],2)` ) and `r round(BTN$hind[1],2)` (SD = `r round(BTN$hindsd[1],2)` ) for the SEBS and NEBS, respectively.

Model projected trends for sea surface temperature (SST) are similar to those of bottom temperature but are warmer and with higher agreement for all three models under low carbon mitigation scenarios (ssp585). Under low carbon mitigation scenarios, estimates of end of century warming of SST [(2080-2100)-(`r paste0(ref_years[1],"-", rev(ref_years)[1])`)] range from `r a <- SST[SST$scen=="ssp126",]; paste0("+",round((a$mn-a$sd)-a$hind,2)," to +",round((a$mn+a$sd)-a$hind,2))` $^o$C for high carbon mitigation (ssp126) scenarios and `r a <- SST[SST$scen=="ssp585",]; paste0("+",round((a$mn-a$sd)-a$hind,2)," to +",round((a$mn+a$sd)-a$hind,2))` $^o$C for low carbon mitigation (ssp585) scenarios. For the NEBS, estimates of end of century warming of SST [(2080-2100)-(`r paste0(ref_years[1],"-", rev(ref_years)[1])`)] range from `r a <- SSTN[SSTN$scen=="ssp126",]; paste0("+",round((a$mn-a$sd)-a$hind,2)," to +",round((a$mn+a$sd)-a$hind,2))` $^o$C for high carbon mitigation (ssp126) scenarios and `r a <- SSTN[SSTN$scen=="ssp585",]; paste0("+",round((a$mn-a$sd)-a$hind,2)," to +",round((a$mn+a$sd)-a$hind,2))` $^o$C for low carbon mitigation (ssp585) scenarios (Fig. 2). Mean historical (`r paste0(ref_years[1],"-", rev(ref_years)[1])`) SSTs from the hindcast are `r round(SST$hind[1],1)` (SD = `r round(SST$hindsd[1],2)` ) and `r round(SSTN$hind[1],2)` (SD = `r round(SSTN$hindsd[1],2)` ) for the SEBS and NEBS, respectively.


Global Warming Levels (GWL) of +3 and +4 $^o$C over pre-industrial temperatures (1850-1900) are associated with significant warming of surface and bottom water temperatures in both the Northern Bering and Southern Bering Seas (NEBS and SEBS, respectively; Fig. 3). In contrast, under GWL of +1.5 $^o$C, SST and BT in both regions remain within the near-present range of climate variability (2010-2021). GWL of +1.5 and +2 $^o$C represent the target and limit respectively of the UNFCCC [Paris Agreement and Nationally determined contributions (NDCs)](https://unfccc.int/ndc-synthesis-report-2022#Projected-GHG-Emission-levels); UNFCCC 2022), i.e. the GWLs beyond which climate change impacts and risks increase,- and feasibility and effectiveness of adaptation actions decrease-, rapidly with each GWL (IPCC 2022). As a point of reference, in the most recent 6th assessment report the IPCC found that 2019 GWL was +1.08 $^o$C (IPCC 2021). 

Both the SEBS and NEBS exhibit seasonal patterns in BT and SST, and seasonally-specific patterns in warming (relative to 1980-2013 climatology from corresponding historical runs for each ESM) across IPCC scenarios and among ESM models (Figs. 4, 5; note these plots are of non-bias corrected values). In general there is agreement in warming trends among the three ensemble members for all seasons in the SEBS, especially under low carbon mitigation (high greenhouse gas emissions; ssp585) scenarios (Fig. 4), although the magnitude varies with model (cesm > miroc > gfdl). In the NEBS there is more variability among the three ensemble members under high mitigation scenarios (ssp126), with 2 of 3 models projecting little warming in winter months(exception is cesm). However, under low carbon mitigation (high greenhouse gas emissions) scenarios (ssp585), all three models project warming in winter months (i.e., reduced sea ice) as well as increases in spring, summer, and fall BT (Fig. 5). 

<!-- [add pH details] -->

## Factors influencing observed trends

For more information about climate change impacts, risks, adaptation, and mitigation see the Intergovernmental Panel on Climate Change (IPCC) [6th Assessment Report (www.ipcc.ch/assessment-report/ar6/)](https://www.ipcc.ch/assessment-report/ar6/) and [www.climate.gov](www.climate.gov).

Carbon dioxide (CO$_2$) is naturally occurring greenhouse gas (GHG) that along with other GHGs acts to absorb and re-emit infrared energy (heat) from solar radiation, warming the earth's surface (i.e., the ‘greenhouse effect’). Naturally occurring CO$_2$ (ocean off-gassing, volcanoes, etc.) is naturally offset by carbon sinks (e.g., photosynthesis of plants on land and in the ocean) acting to keep CO$_2$ relatively stable for more than 800,000 years at or below 300 parts per million (ppm). However, scientific observations and models have shown that atmospheric CO$_2$ concentrations have been rising steadily over the past century due to anthropogenic (human) activities, primarily the the burning of fossil fuels for energy and other uses. Rates of anthropogenic CO$_2$ release into the atmosphere exceed natural carbon sinks and have resulted in a rapid accumulation of atmospheric CO$_2$. The IPCC 6th Assessment Report states that *"observed increases in well-mixed greenhouse gas (GHG) concentrations since around 1750 are unequivocally caused by human activities"* and that the *"land and ocean have taken up a near-constant proportion (globally about 56% per year) of CO$_2$ emissions from human activities over the past six decades, with regional differences (high confidence)"* (IPCC 2021).

Current atmospheric CO$_2$ levels of 410 ppm (IPCC 2021) have not been experienced for at least 2 million years, and the rate of increase in CO$_2$ over the last century is unprecedented in the last 800,000 years (based on multiple lines of evidence including Antarctic ice core data and isotopes IPCC 2021). There is a near-linear relationship between cumulative CO$_2$ emissions and increases in global surface temperature (IPCC 2021), and changes in atmospheric CO$_2$ and associated warming have direct impacts on ocean processes and chemistry. In the most recent assessment report, the IPPC (2021) states *"better integration of paleo-oceanographic data with modelling along with higher-resolution analyses of transient changes have improved understanding of long-term ocean processes...This paleo context supports the assessment that ongoing increase in ocean heat content (OHC) represents a long-term commitment, essentially irreversible on human time scales (high confidence)"*. This paleo context has also helped illuminate the complex role of oceans in the regulation of the global climate and atmospheric CO$_2$ during previous glacial–interglacial warming intervals.  Presently, absorption of atmospheric heat by the world's oceans increases ocean temperatures, warming surface waters to depth and raising the "ocean heat content". Absorption of atmospheric CO$_2$ alters the chemistry of the ocean increasing acidity and lowering the pH. In addition, atmospheric warming alters physical and chemical processes (e.g., precipitation, wind patterns, sea level, ocean circulation, and sea ice thickness and extent) in ways that further change the ocean and atmospheric cycles, i.e., the climate of a given region. Accordingly, the IPCC (2021) states, "it is *virtually certain* that the global upper ocean (0–700 m) has warmed since the 1970s and *extremely likely* that human influence is the main driver. It is *virtually certain* that human-caused CO$_2$ emissions are the main driver of current global acidification of the surface open ocean". 

The near linear relationship between cumulative CO$_2$ emissions and increases in global surface temperature (IPCC 2021) has enabled scientists evaluate future climate conditions under alternative CO$_2$ and GHG emissions scenarios, known as Shared Socioeconomic Pathways (SSPs; O'neil et al. 2017). These allow for projections of changes in climate and ocean temperature and chemistry under Global Warming Levels (GWLs). Patterns of warming reported in this contribution reflect global changes in atmospheric carbon, climate conditions, and oceanic conditions from Earth Systems Models, but are refined through a regional lens via a the high resolution Bering10K ROMSNPZ ocean model that is able to replicate fine scale oceanographic processes (e.g., changes in sea ice and circulation on a short timestep) that act to amplify or attenuate larger scale climate change effects.

## Implications: 

Historically, warming temperatures and marine heat waves have been associated with changes to food-web dynamics, species redistribution, and ecosystem structure and processes (Huntington et al. 2020). Projected ocean warming from global models is associated with declines in marine fish biomass, benthic biomass, and fisheries catch potential (IPCC AR6 WGII). Evaluations of projected temperature effects in Bering sea ecosystems and fisheries under high emission scenarios (ssp585) are still underway but include modeled declines in winter sea ice and summer cold pool extent associated with increased warming under low carbon mitigation scenarios (ssp585). Increased warming in EBS projections is also associated with emergent declines in modeled fall euphausiid and large copepod biomass (Hermann et al. 2021), shifts in spring bloom timing to earlier (30-60 d) and slightly larger phytoplankton and zooplankton blooms (relative to hindcasts), and declines in the magnitude of fall total phytoplankton and large zooplankton blooms (Cheng et al. 2021).

## Description of index: 

We report trends in modeled bottom temperature and sea surface temperature from a 30-layer Bering Sea regional oceanographic model at 10km horizontal resolution which has incorporated lower trophic level biology (Kearney et al. 2020) and marine carbonate chemistry (Pilcher et al. 2019). See the [*Bering 10K dataset documentation*](https://zenodo.org/record/4586950/files/Bering10K_dataset_documentation.pdf) for more information and technical details. We present the [*Alaska NOAA Integrated Ecosystem Assessment (IEA) program*](www.integratedecosystemassessment.noaa.gov) annual hindcast and two CMIP6 emission scenarios projected as part of the [*Alaska Climate Integrated Modeling (ACLIM) project*](www.fisheries.noaa.gov/alaska/ecosystems/alaska-climate-integrated-modeling-project). In this, a low atmospheric carbon emissions scenario (ssp126) and a low carbon mitigation scenario (ssp585; O’Neil et al. [2017](
https://link.springer.com/article/10.1007/s10584-013-0905-2)) and three global Earth System Models (ESMs; 'cesm', 'gfdl' and 'mir') were selected from the [Coupled Model Intercomparison Project (CMIP6)](https://www.wcrp-climate.org/wgcm-cmip/wgcm-cmip6) and used to force the regional model. See Hermann et al. 2021, Cheng et al. 2021, Kearney et al. 2020, and Pilcher et al. 2019 for details about the regional model projections and Hollowed et al. 2020 for details about the [ACLIM](www.fisheries.noaa.gov/alaska/ecosystems/alaska-climate-integrated-modeling-project) project and forcing (climate scenario and ESM) selection.

In support of [ACLIM](www.fisheries.noaa.gov/alaska/ecosystems/alaska-climate-integrated-modeling-project), a number of different biophysical index timeseries were calculated based on the Bering10K simulations and provide the primary means of linking the physical and lower trophic level dynamics to the ACLIM suite of upper trophic level and socioeconomic models; see Hollowed et al. 2020 for further details. The timeseries reported here are derived from the area-weighted strata averages for Summer (months Jul - Sep) and Winter (Jan - Feb) for the Northern Bering Sea (strata `r paste(NEBS_strata,collapse=", ")` of the eastern Bering Sea shelf bottom trawl survey of the Alaska Fisheries Science Center) and Southern Bering sea (strata `r  paste(SEBS_strata,collapse=", ")`). The timerseries were bias-corrected to hindcast simulations using historical forcing (during 1980-2013) from each ESM. More detail on this approach is available by request.

The climate simulations presented here are dynamically downscaled from a selection of the historical and shared socioeconomic pathway simulations from the sixth phase of the [CMIP6](https://www.wcrp-climate.org/wgcm-cmip/wgcm-cmip6). Names reflect the parent global model simulation (miroc = MIROC ES2L, cesm = CESM2, gfdl = GFDL ESM4) and emissions scenario via Shared Socioeconomic Pathways (SSPs) (ssp126 = SSP1-2.6, ssp585 = SSP5-8.5, historical = Historical). ssp126 represents a high carbon mitigation (low greenhouse gas emissions) scenario; ssp585 represents the low carbon mitigation scenario. More information on the SSPs and their use in climate projections is available at O’Neil et al. [2017](
https://link.springer.com/article/10.1007/s10584-013-0905-2).

To determine mean temperatures associated with standardized levels of global warming for each scenario and ESM, we used CMIP6 Global Warming Levels from Hauser et al. [2021](https://doi.org/10.5281/zenodo.4600706) and publicly available at [https://github.com/mathause/cmip_warming_levels](https://github.com/mathause/cmip_warming_levels).


## Literature Cited

Cheng et al., 2021 W. Cheng, A.J. Hermann, A.B. Hollowed, K.K. Holsman, K.A. Kearney, D.J. Pilcher, et al. Eastern Bering Sea shelf environmental and lower trophic level responses to climate forcing: results of dynamical downscaling from CMIP6 Deep-Sea Res. II, 193 (2021), Article 104975, [10.1016/j.dsr2.2021.104975](https://www.sciencedirect.com/science/article/pii/S0967064521000515).

Hauser, M., F. Engelbrecht, F. & E. Fischer (2021). Transient global warming levels for CMIP5 and CMIP6 (v0.2.0) [Data set]. https://github.com/mathause/cmip_warming_levels/tree/v0.2.0. [doi.org/10.5281/zenodo.4600706](https://doi.org/10.5281/zenodo.4600706).

Hermann et al., 2021 A.J. Hermann, K. Kearney, W. Cheng, D. Pilcher, K. Aydin, K.K. Holsman, et al. Coupled modes of projected regional change in the Bering Sea from a dynamically downscaling model under CMIP6 forcing Deep-Sea Res. II (2021), [10.1016/j.dsr2.2021.104974 194 104974](https://www.sciencedirect.com/science/article/pii/S0967064521000503).

Hollowed, K. K. Holsman, A. C. Haynie, A. J. Hermann, A. E. Punt, K. Y. Aydin, J. N. Ianelli, S. Kasperski, W. Cheng, A. Faig, K. Kearney, J. C. P. Reum, P. D. Spencer, I. Spies, W. J. Stockhausen, C. S. Szuwalski, G. Whitehouse, and T. K. Wilderbuer. Integrated modeling to evaluate climate change impacts on coupled social-ecological systems in Alaska. Frontiers in Marine Science, 6(January):1–18, 2020. DOI: [10.3389/fmars.2019.00775](https://www.frontiersin.org/articles/10.3389/fmars.2019.00775/full).

IPCC, 2021: Summary for Policymakers. In: Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change [Masson-Delmotte, V., P. Zhai, A. Pirani, S.L. Connors, C. Péan, S. Berger, N. Caud, Y. Chen, L. Goldfarb, M.I. Gomis, M. Huang, K. Leitzell, E. Lonnoy, J.B.R. Matthews, T.K. Maycock, T. Waterfield, O. Yelekçi, R. Yu, and B. Zhou (eds.)]. Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA, pp. 3−32, [doi:10.1017/9781009157896.001](https://www.ipcc.ch/report/ar6/wg1/downloads/report/IPCC_AR6_WGI_SPM.pdf).

IPCC, 2022: Climate Change 2022: Impacts, Adaptation and Vulnerability. Contribution of Working Group II to the Sixth Assessment Report of the
Intergovernmental Panel on Climate Change [H.-O. Pörtner, D.C. Roberts, M. Tignor, E.S. Poloczanska, K. Mintenbeck, A. Alegría, M. Craig, S. Langsdorf,
S. Löschke, V. Möller, A. Okem, B. Rama (eds.)]. Cambridge University Press. Cambridge University Press, Cambridge, UK and New York, NY, USA,
3056 pp., [doi:10.1017/9781009325844](https://report.ipcc.ch/ar6/wg2/IPCC_AR6_WGII_FullReport.pdf).

Kearney, K,  A. Hermann, W. Cheng, I. Ortiz, and K. Aydin. A coupled pelagic benthic-sympagic biogeochemical model for the Bering Sea: documentation and validation of the BESTNPZ model (v2019.08.23) within a high resolution regional ocean model. Geoscientific Model Development, 13 (2):597–650, 2020. DOI: [10.5194/gmd13-597-2020](https://www.geosci-model-dev.net/13/597/2020/10 https://github.com/beringnpz/roms-bering-sea).

O’Neill, B. C.,E. Kriegler, K. Riahi, K. L. Ebi, S. Hallegatte, T. R. Carter, R. Mathur & D. P. van Vuuren (2017) A new scenario framework for climate change research: the concept of shared socioeconomic pathways. Climatic Change (2014) 122:387–400 [DOI 10.1007/s10584-013-0905-2](https://link.springer.com/article/10.1007/s10584-013-0905-2)

Pilcher, D. J., D. M. Naiman, J. N. Cross, A. J. Hermann, S. A. Siedlecki, G. A. Gibson, and J. T. Mathis. Modeled Effect of Coastal Biogeochemical Processes, Climate Variability, and Ocean Acidification on Aragonite Saturation State in the Bering Sea. Frontiers in Marine Science, 5(January):1–18, 2019. DOI: [10.3389/fmars.2018.00508](https://www.frontiersin.org/articles/10.3389/fmars.2018.00508/full) 12.

Pilcher, D. J.,  J.N.Cross, A.J.Hermanna, K.A.Kearneya, W.Chenga, J.T.Mathis. 2022. Dynamically downscaled projections of ocean acidification for the Bering Sea. Deep Sea Res II. (198) [https://doi.org/10.1016/j.dsr2.2022.105055](https://doi.org/10.1016/j.dsr2.2022.105055).

UNFCCC 2022. Nationally determined contributions under the Paris Agreement Synthesis report by the secretariat. United Nations Framework Convention on Climate Change. Conference of the Parties serving as the meeting of the Parties to the Paris Agreement Fourth session. Sharm el-Sheikh, 6–18 November 2022. FCCC/PA/CMA/2022/4 [https://unfccc.int/ndc-synthesis-report-2022](https://unfccc.int/ndc-synthesis-report-2022).

## Figures: 

![Figure 1. Bias corrected summer sea surface  temperature (top row) and bottom temperature (bottom row) for the southern Bering Sea (SEBS) from the hindcast (dark blue line) and projections under low (ssp126, left column; cool colors) and high (ssp585, right column, warm colors) emission scenarios; individual Earth System Models are shown as individual lines. Average modeled temperatures from the reference period (1980 -2013) of the hindcast are show as the horizontal blue line; dashed lines represent +/- 1 standard deviation of the mean. Note different scales between rows.](ESR_EBS/Figs/ACLIM_Holsman_Fig1.png){ width=100% }

![Figure 2, Bias corrected summer sea surface  temperature (top row) and bottom temperature (bottom row) for the Northern Bering Sea (NEBS) from the hindcast (dark blue line) and projections under low (ssp126, left column; cool colors) and high (ssp585, right column, warm colors) emission scenarios; individual Earth System Models are shown as individual lines. Average modeled temperatures from the reference period (1980 -2013) of the hindcast are show as the horizontal blue line; dashed lines represent +/- 1 standard deviation of the mean. Note different scales between rows.](ESR_EBS/Figs/ACLIM_Holsman_Fig2.png){ width=100% }

![Figure 3. Southern and Northern Bering Sea ('SEBS' and 'NEBS', respectively) modeled summer bottom and sea surface temperatures ('BT' and 'SST', respectively) as a function of CMIP6 Global Warming Levels (mean global increase in temperature relative to pre-industrial temperatures (1850-1900)). Recent hindcast ranges are reported ('2010-2021') as well as bias corrected projections from the Bering10K model for each GWL (+1 to +4 GWL). Boxplots represent the 25th and 75th percentile (i.e, the interquartile range) with the horizontal line representing the median temperature, and the error bars representing the min or max (IQR +/- IQR*1.5). Outliers are represented by points (e.g., MHW years if above the boxplot). For more information on interpretation of boxplots see https://r-graph-gallery.com/boxplot.html](ESR_EBS/Figs/ACLIM_Holsman_Fig3.png){ width=100% }

![Figure 4. Southern Bering Sea (SEBS) bottom water temperature (oC) projected under two climate scenarios; high carbon mitigation via Shared Socioeconomic Pathways (ssp126, left column); low carbon mitigation, (ssp585, right column). Rows reflect the parent global model simulation (miroc = MIROC ES2L, cesm = CESM2, gfdl = GFDL ESM4) dynamically downscaled to a high resolution regional model (Bering10K K20P19 30 layer ROMSNPZ model). Average modeled climatology from the reference period (1980 -2013) of the historical run for each ESM is show as the solid blue line; dashed lines represent +/- 1 standard deviation of the mean.](ESR_EBS/Figs/ACLIM_Holsman_Fig4.png){ width=100% }

![Figure 5. Northern Bering Sea (NEBS) bottom water temperature (oC) projected under two climate scenarios; high carbon mitigation via Shared Socioeconomic Pathways (ssp126, left column); low carbon mitigation, (ssp585, right column). Rows reflect the parent global model simulation (miroc = MIROC ES2L, cesm = CESM2, gfdl = GFDL ESM4) dynamically downscaled to a high resolution regional model (Bering10K K20P19 30 layer ROMSNPZ model). Average modeled climatology from the reference period (1980 -2013) of the historical run for each ESM is show as the solid blue line; dashed lines represent +/- 1 standard deviation of the mean.](ESR_EBS/Figs/ACLIM_Holsman_Fig5.png){ width=100% }

<!-- ![Bottom water pH pojected under two climate scenarios ( high carbon mitigation (ssp126); low carbon mitigation, (ssp585) from three General Circulation Models (GCMs) dynamically downscaled to a high resolution regional model (Bering10K 30 layer ROMSNPZ model).](ESR_EBS/Figs/nonBC_weeklyProj_pH.png){ width=100% } -->

<!-- ![Bottom water pH pojected under two climate scenarios ( high carbon mitigation (ssp126); low carbon mitigation, (ssp585) from three General Circulation Models (GCMs) dynamically downscaled to a high resolution regional model (Bering10K 30 layer ROMSNPZ model).](ESR_EBS/Figs/nonBC_weeklyProj_pH_N.png){ width=100% } -->