Merge pull request #33 from kdorheim/fix_ohc

Try optimizing for OHC

results/12_alpha_ecs_ohc_unc_nmse.txt

@@ -0,0 +1,47 @@
+parameters	values
+beta	0.64953154546281
+q10_rh	1.76
+diff	1.042
+S	2.32537207553702
+alpha	0.438414955305681
+
+Objective Function Value: 4.77e-05
+Counts: 40
+Counts: 40
+Convergence:  0
+Messages: CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH
+
+CO2 MSE: 8.57
+T MSE:   0.0127
+RMSE:    2.93
+T MSE accounting for unc: 0.0026
+OHC MSE accounting for unc: 9.71
+
+***Key Metrics***
+TCRE: 1.23
+TCR:  1.67
+
+***Historical Warming and ERF***
+GSAT Warming:              1.01
+Ocean Heat Content Change: 425
+Total Aerosol ERF:         -1.24
+WMGHG ERF:                 3.78
+Methane ERF:               0.545
+
+***Future Warming***
+scenario	start	end	GSAT
+ssp119	2021	2040	0.51
+ssp119	2041	2060	0.484
+ssp119	2081	2100	0.138
+ssp126	2021	2040	0.545
+ssp126	2041	2060	0.688
+ssp126	2081	2100	0.473
+ssp245	2021	2040	0.603
+ssp245	2041	2060	1
+ssp245	2081	2100	1.4
+ssp370	2021	2040	0.654
+ssp370	2041	2060	1.24
+ssp370	2081	2100	2.51
+ssp585	2021	2040	0.7
+ssp585	2041	2060	1.43
+ssp585	2081	2100	3.16

results/default_stats.txt

@@ -6,6 +6,15 @@ Unnamed Hector core	NA	diff	2.38	cm2/s
 Mean of T, CO2 MSEs: 2.31
 Mean of T, CO2 NMSEs: 0.000948
 Mean of T, CO2 NMSEs (with unc): 0.000313
+scenario	year	variable	value	units
+Unnamed Hector core	NA	beta	0.53	(unitless)
+Unnamed Hector core	NA	q10_rh	1.76	(unitless)
+Unnamed Hector core	NA	diff	2.38	cm2/s
+
+Mean of T, CO2 MSEs: 2.31
+Mean of T, CO2 NMSEs: 0.000948
+Mean of T, CO2 NMSEs (with unc): 0.000313
+Mean of T, CO2, OHC NMSEs (with unc): 0.00022
 
 Mean of smooth T, CO2 MSEs (k = 3 ): 2.3
 Mean of smooth T, CO2 NMSEs (k = 3 ): 0.000864
@@ -22,6 +31,7 @@ RMSE:    2.14
 CO2 NMSE: 2.74e-07
 T NMSE:   0.0019
 T NMSE with unc: 0.000625
+OHC NMSE with unc: 3.43e-05
 
 ***Key Metrics***
 TCRE: 1.51

scripts/alpha_ohc_unc_nmse.R

@@ -0,0 +1,88 @@
+# Script for experiment 12
+# Script to use normalized T, OHC and CO2 MSEs while accounting for T, OHC unc
+# Also includes ECS and alpha as params to optimize over
+# Uses OHC range from Matilda table rather than manuscript script
+# Author: Peter Scully
+# Date: 6/25/24
+
+### Constants and Imports ###
+
+# Importing libraries
+library(hector)
+
+# Setting up file paths
+COMP_DATA_DIR <- file.path(here::here(), "comparison_data")
+SCRIPTS_DIR <- file.path(here::here(), "scripts")
+RESULTS_DIR <- file.path(here::here(), "results")
+
+CO2_PATH <- file.path(COMP_DATA_DIR,
+                      "Supplementary_Table_UoM_GHGConcentrations-1-1-0_annualmeans_v23March2017.csv")
+TEMP_PATH <-
+  file.path(COMP_DATA_DIR,
+            "HadCRUT.5.0.2.0.analysis.summary_series.global.annual.csv")
+
+OHC_PATH <- file.path(COMP_DATA_DIR, "OHC_ensemble_Kuhlbrodt_etal_2022.csv")
+
+
+INI_FILE <- system.file("input/hector_ssp245.ini", package = "hector")
+PARAMS <- c(BETA(), Q10_RH(), DIFFUSIVITY(), ECS(), AERO_SCALE())
+
+OUTPUT <- file.path(RESULTS_DIR, "12_alpha_ecs_ohc_unc_nmse.txt")
+
+
+source(file.path(SCRIPTS_DIR, "major_functions.R"))
+
+### Getting observational data ###
+co2_data <- get_co2_data(CO2_PATH, include_unc = TRUE)
+temp_data <- get_temp_data(TEMP_PATH, include_unc = TRUE)
+ohc_data <- get_ohc_data(OHC_PATH, include_unc = T)
+obs_data <- rbind(co2_data, temp_data, ohc_data)
+
+### Calling optim ###
+best_pars <- run_optim(obs_data = obs_data,
+                       ini_file = INI_FILE,
+                       params = PARAMS,
+                       lower = c(0.5 - 0.232, 2.2 - 0.44, 1.16 - 0.118, 2, 0),
+                       upper = c(0.5 + 0.232, 2.2 + 0.44, 1.16 + 0.118, 5, 3),
+                       yrs = 1750:2014,
+                       vars = c(GMST(), CONCENTRATIONS_CO2(), HEAT_FLUX()),
+                       error_fn = mean_T_CO2_OHC_nmse_unc,
+                       include_unc = T,
+                       method = "L-BFGS-B",
+                       output_file = OUTPUT)
+
+### Outputting individual MSEs ###
+hector_data <- run_hector(ini_file = INI_FILE, 
+                          params = PARAMS, 
+                          vals = best_pars, 
+                          yrs = 1750:2014, 
+                          vars = c(GMST(), CONCENTRATIONS_CO2(), HEAT_FLUX()))
+
+T_mse <- get_var_mse(obs_data = obs_data, 
+                     hector_data = hector_data, 
+                     var = GMST(), 
+                     yrs = 1850:2014)
+CO2_mse <- get_var_mse(obs_data = obs_data, 
+                       hector_data = hector_data, 
+                       var = CONCENTRATIONS_CO2(), 
+                       yrs = c(1750, 1850:2014))
+T_mse_unc <- get_var_mse_unc(obs_data = obs_data,
+                             hector_data = hector_data,
+                             var = GMST(),
+                             yrs = 1850:2014,
+                             mse_fn = mse_unc)
+OHC_mse_unc <- get_var_mse_unc(obs_data = obs_data, 
+                                hector_data = hector_data, 
+                                var = "OHC", 
+                                yrs = 1957:2014,
+                                mse_fn = mse_unc)
+
+write_metric("CO2 MSE:", CO2_mse, OUTPUT)
+write_metric("T MSE:  ", T_mse, OUTPUT)
+write_metric("RMSE:   ", sqrt(mean(CO2_mse, T_mse)), OUTPUT) # not 100% sure this is how we want to calculate this
+write_metric("T MSE accounting for unc:", T_mse_unc, OUTPUT)
+write_metric("OHC MSE accounting for unc:", OHC_mse_unc, OUTPUT)
+write("", OUTPUT, append = TRUE)
+
+### Outputting table metrics ###
+calc_table_metrics(PARAMS, best_pars, OUTPUT)

scripts/default_stats.R

@@ -19,6 +19,8 @@ TEMP_PATH <-
   file.path(COMP_DATA_DIR,
             "HadCRUT.5.0.2.0.analysis.summary_series.global.annual.csv")
 
+OHC_PATH <- file.path(COMP_DATA_DIR, "OHC_ensemble_Kuhlbrodt_etal_2022.csv")
+
 INI_FILE <- system.file("input/hector_ssp245.ini", package = "hector")
 PARAMS <- c(BETA(), Q10_RH(), DIFFUSIVITY())
 
@@ -31,7 +33,8 @@ source(file.path(SCRIPTS_DIR, "major_functions.R"))
 ### Getting observational data ###
 co2_data <- get_co2_data(CO2_PATH, include_unc = T)
 temp_data <- get_temp_data(TEMP_PATH, include_unc = T)
-obs_data <- rbind(co2_data, temp_data)
+ohc_data <- get_ohc_data(OHC_PATH, include_unc = T)
+obs_data <- rbind(co2_data, temp_data, ohc_data)
 
 
 ### Outputting default values ###
@@ -53,7 +56,7 @@ hector_data <- run_hector(ini_file = INI_FILE,
                           params = NULL, 
                           vals = best_pars, 
                           yrs = 1750:2014, 
-                          vars = c(GMST(), CONCENTRATIONS_CO2()),
+                          vars = c(GMST(), CONCENTRATIONS_CO2(), HEAT_FLUX()),
                           include_unc = T)
 
 
@@ -67,6 +70,9 @@ write_metric("Mean of T, CO2 NMSEs:",
 write_metric("Mean of T, CO2 NMSEs (with unc):",
              mean_T_CO2_nmse_unc(obs_data, hector_data),
              OUTPUT)
+write_metric("Mean of T, CO2, OHC NMSEs (with unc):",
+             mean_T_CO2_OHC_nmse_unc(obs_data, hector_data),
+             OUTPUT)
 write("", OUTPUT, append = TRUE)
 
 # Outputting smoothed MSEs
@@ -122,6 +128,11 @@ CO2_nmse <- get_var_mse(obs_data = obs_data,
                        var = CONCENTRATIONS_CO2(), 
                        yrs = c(1750, 1850:2014),
                        mse_fn = nmse)
+OHC_nmse_unc <- get_var_mse_unc(obs_data = obs_data, 
+                                hector_data = hector_data, 
+                                var = "OHC", 
+                                yrs = 1957:2014,
+                                mse_fn = nmse_unc)
 
 write_metric("CO2 MSE:       ", CO2_mse, OUTPUT)
 write_metric("T MSE:         ", T_mse, OUTPUT)
@@ -131,6 +142,7 @@ write("", OUTPUT, append = TRUE)
 write_metric("CO2 NMSE:", CO2_nmse, OUTPUT)
 write_metric("T NMSE:  ", T_nmse, OUTPUT)
 write_metric("T NMSE with unc:", T_nmse_unc, OUTPUT)
+write_metric("OHC NMSE with unc:", OHC_nmse_unc, OUTPUT)
 write("", OUTPUT, append = TRUE)
 
 ### Outputting table metrics ###

scripts/error_functions.R

@@ -269,6 +269,34 @@ mean_T_CO2_nmse_unc <- function(obs_data, hector_data) {
   return(mean(c(T_mse, CO2_mse)))
 }
 
+# mean_T_CO2_OHC_nmse_unc: function to find the mean of the temperature, CO2, 
+#                          and OHC NMSEs between observed and predicted data
+#                          while accounting for temperature and OHC uncertainty
+#
+# args: 
+#   obs_data    - data frame of observed data formatted like Hector data frame
+#   hector_data - data frame outputted by Hector
+#
+# Returns: Average NMSE between predicted and observed data
+mean_T_CO2_OHC_nmse_unc <- function(obs_data, hector_data) {
+  T_mse <- get_var_mse_unc(obs_data = obs_data, 
+                           hector_data = hector_data, 
+                           var = GMST(), 
+                           yrs = 1850:2014,
+                           mse_fn = nmse_unc)
+  CO2_mse <- get_var_mse(obs_data = obs_data, 
+                         hector_data = hector_data, 
+                         var = CONCENTRATIONS_CO2(), 
+                         yrs = c(1750, 1850:2014),
+                         mse_fn = nmse)  
+  OHC_mse <- get_var_mse_unc(obs_data = obs_data,
+                             hector_data = hector_data,
+                             var = "OHC",
+                             yrs = 1957:2014,
+                             mse_fn = nmse_unc)
+  return(mean(c(T_mse, CO2_mse, OHC_mse)))
+}
+
 # smooth_T_CO2_mse: function to find the mean of smoothed temperature and CO2
 #                   MSEs between observed & predicted data for a given variable
 #

scripts/graph_comparison_plots.R

@@ -88,8 +88,15 @@ alpha_bb_data <- run_hector(ini_file = INI_FILE,
                             vars = c(GMST(), CONCENTRATIONS_CO2()))
 alpha_bb_data$scenario <- "Hector - NMSE w/ unc, Big Box Fit w/ S, Alpha"
 
+alpha_ohc_data <- run_hector(ini_file = INI_FILE,
+                         params = PARAMS,
+                         vals = c(0.65, 1.76, 1.04, 2.33, 0.438),
+                         yrs = 1750:2014,
+                         vars = c(GMST(), CONCENTRATIONS_CO2()))
+alpha_ohc_data$scenario <- "Hector - NMSE w/ unc (incl. OHC), Fit w/ S, Alpha"
+
 #hector_data <- rbind(default_data, nmse_data, nmse_bb_data, ecs_data, ecs_bb_data, alpha_data, alpha_bb_data)
-hector_data <- rbind(default_data, alpha_data)
+hector_data <- rbind(default_data, alpha_data, alpha_ohc_data)
 hector_data$lower <- hector_data$value
 hector_data$upper <- hector_data$value
 

scripts/major_functions.R

@@ -98,6 +98,49 @@ get_temp_data <- function(file, scenario = "historical", include_unc = F) {
   return(temp_data)
 }
 
+# get_ohc_data - function to get ocean heat content data
+#
+# args:
+#   file - path to historical OHC data file
+#   scenario - name of scenario being run (default: "historical")
+#   include_unc - boolean indicating whether to include uncertainty data
+#                 (default: F)
+#
+# returns: Hector-style data frame with OHC data
+get_ohc_data <- function(file, scenario = "historical", include_unc = F) {
+
+  # Reading in only OHC data
+  ohc_data <- read.table(file, 
+                         skip = 2, 
+                         sep = ",",
+                         colClasses = c("numeric", "NULL", "NULL", "NULL",
+                                        "NULL", "NULL", "NULL", "numeric",
+                                        "numeric"))
+
+  # Fixing table formatting
+  ohc_data <- na.omit(ohc_data)
+  colnames(ohc_data) <- c("year", "value", "unc")
+
+  # Getting rid of non-integer years
+  ohc_data$year <- ohc_data$year - 0.5
+
+  # Adding in new columns to match Hector data frames
+  ohc_data$scenario <- scenario
+  ohc_data$variable <- "OHC"
+  ohc_data$units <- "ZJ"
+
+  # Adding in confidence interval (if applicable)
+  if (include_unc) {
+    ohc_data$lower <- ohc_data$value - ohc_data$unc
+    ohc_data$upper <- ohc_data$value + ohc_data$unc
+  }
+
+  # Getting rid of raw uncertainty column
+  ohc_data$unc <- NULL
+
+  return(ohc_data)
+}
+
 
 ############################
 ### Data Frame Functions ###
@@ -257,6 +300,31 @@ run_hector <- function(ini_file, params, vals, yrs, vars, include_unc = F) {
     data <- rel_to_interval(data = data, var = GMST(), start = 1961, end = 1990)
   }
 
+  # Measuring ocean heat content (if applicable)
+  if (HEAT_FLUX() %in% vars) {
+
+    ### Converting heat fluxes to cumulative heat content
+
+    # Getting heat flux data
+    ohc_data <- filter(data, variable == HEAT_FLUX() & year %in% 1957:2014)
+
+    # Converting heat flux to OHC change by year
+    ohc_data$value <- ohc_data$value * OCEAN_AREA * W_TO_ZJ
+
+    # Converting OHC changes by year to total OHC (relative to 1957)
+    ohc_data$value <- cumsum(ohc_data$value)
+    ohc_data$variable <- "OHC"
+
+    # Making OHC relative to 2005-2014 average
+    ohc_data <- rel_to_interval(data = ohc_data,
+                                 var = "OHC",
+                                 start = 2005,
+                                 end = 2014)
+
+    # Adding ohc data to Hector data frame
+    data <- rbind(data, ohc_data)
+  }
+
   # Adding in upper and lower bounds (if applicable)
   if (include_unc) {
     data$upper <- data$value