Hector: GMST vs GMAT

Objective

RCMIP II is interested in the what they refer to as the Surface Air Ocean Blended Temperature Change aka in other words the GMST instead of the GMAT. Currently (September 2020) Hector does not return that as an output but it is possible to calculate from different output streams. Which is what we do here in preparation for the format netcdf steps. The objective of this document is to back calculate the GMST from Hector.

Reasons why this might matter

• we may want to opt for using GMST in calibration process & impact the calibration protocol
  
• may want to add it as standard output

Set Up and Run Hector

library(dplyr)
library(tidyr)
library(ggplot2)
library(data.table)
library(hector)
library(magrittr)

Run Hector with the rcp 45 scenario.

# Define ini file / scenario to run. 
ini  <- system.file(package = 'hector', 'input/hector_rcp45.ini')
core <- newcore(ini)
run(core)

Hector core:    unnamed hector core
Start date: 1745
End date:   2300
Current date:   2300
Input file: /Library/Frameworks/R.framework/Versions/3.6/Resources/library/hector/input/hector_rcp45.ini

# Get all of the different temperature related outputs. 
output <- fetchvars(core, 1750:2100, vars = c(GLOBAL_TEMP(), LAND_AIR_TEMP(), OCEAN_AIR_TEMP()))
output <- as.data.table(output)

Calculate GMST

We know based off of a conversation with Steve

Right now hector reports global temperature as: temp[tstep] = flnd * temp_landair[tstep] + (1.0 - flnd) * bsi * temp_sst[tstep]; temp_oceanair = bsi * temp_sst[tstep]; So its reporting now surface air temperature as Tgav, but you can report blended temperature by just calculating Tblend = flnd * temp_landair[tstep] + (1.0 - flnd) * temp_sst[tstep]; (since flnd is a physical constant and the other two variables are reported out).

Based on default Hector constants we know that.

flnd = 0.29; 
bsi = 1.3

Because the temp_sst is not read out as an output it has to be back calculated from the ocean air.

# Extract the two different types of air temps 
temp_landair  <- output[variable == LAND_AIR_TEMP(), value]
temp_oceanair <- output[variable == OCEAN_AIR_TEMP(), value]
# Convert the ocean air temp back to sst. 
temp_sst <- temp_oceanair / bsi

#  Calculate the GMST based on the equation that Steve shared. 
GMST_values <- flnd * temp_landair + (1 - flnd) * temp_sst
GMST <- output[variable == LAND_AIR_TEMP(), ] # Copy over the data frame strucutre 
GMST$value <- GMST_values # replace the values 
GMST$variable <- 'GMST'

rslt  <- rbind(GMST, output)
rslt$variable <- ifelse(rslt$variable  == 'Tgav', 'GMAT', rslt$variable)

Visualize Results

ggplot(output) + 
  geom_line(aes(year, value, color  = variable)) + 
  theme_bw()  + 
  labs(y = expression(degree~'C'), 
       title = 'Comparison of Hector Air Temperature Outputs')

As expected the temperature over the land is much warmer than the air temperature over the ocean.

Now let’s compare the global mean air and surface temperature 

ggplot(data = rslt[variable %in% c('GMST', 'GMAT'), ]) + 
  geom_line(aes(year, value, color = variable)) + 
  theme_bw()  + 
  labs(y = expression(degree~'C'), 
       title = 'Comparison of Hector GMST and GMAT')

So now if we look at difference in temperature over time.

GMST_dt <- rslt[variable == 'GMST', ]
GMAT_dt <- rslt[variable == 'GMAT', ]
difference_df <- GMAT_dt
difference_df$variable <- 'GMAT - GMST'
difference_df$value <- GMAT_dt$value - GMST_dt$value

ggplot(data = difference_df) + 
  geom_line(aes(year, value, color = variable)) + 
    theme_bw()  + 
  labs(y = expression(degree~'C'), 
       title = 'GMAT - GMST')

Conclusion

This I think looks good… I think I am somewhat surprised by how similar the historical temperature trends are.

---

Reproducibility

sessionInfo()

R version 3.6.3 (2020-02-29)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Mojave 10.14.6

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] magrittr_1.5      data.table_1.12.8 ggplot2_3.3.2     tidyr_1.1.0       dplyr_1.0.0       hector_2.5.0     

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.4.6     pillar_1.4.4     compiler_3.6.3   base64enc_0.1-3  tools_3.6.3      digest_0.6.25   
 [7] jsonlite_1.6.1   evaluate_0.14    lifecycle_0.2.0  tibble_3.0.1     gtable_0.3.0     pkgconfig_2.0.3 
[13] rlang_0.4.6      rstudioapi_0.11  yaml_2.2.1       xfun_0.14        withr_2.2.0      stringr_1.4.0   
[19] knitr_1.28       generics_0.0.2   vctrs_0.3.1      grid_3.6.3       tidyselect_1.1.0 glue_1.4.1      
[25] R6_2.4.1         rmarkdown_2.3.1  farver_2.0.3     purrr_0.3.4      scales_1.1.1     ellipsis_0.3.1  
[31] htmltools_0.5.0  rsconnect_0.8.16 colorspace_1.4-1 labeling_0.3     stringi_1.4.6    munsell_0.5.0   
[37] crayon_1.3.4    

# Ideally there would be something to figure out which commit of Hector working with