Skip to content

Commit

Permalink
Tried optimizing for OHC
Browse files Browse the repository at this point in the history
  • Loading branch information
ptrscll committed Jun 25, 2024
1 parent ca13527 commit be355ee
Show file tree
Hide file tree
Showing 6 changed files with 146 additions and 2 deletions.
47 changes: 47 additions & 0 deletions results/12_alpha_ecs_ohc_unc_nmse.txt
Original file line number Diff line number Diff line change
@@ -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
Binary file modified results/alpha_comparison_plots.jpeg
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file modified results/ohc_plot.jpeg
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
88 changes: 88 additions & 0 deletions scripts/alpha_ohc_unc_nmse.R
Original file line number Diff line number Diff line change
@@ -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)
9 changes: 8 additions & 1 deletion scripts/graph_comparison_plots.R
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
4 changes: 3 additions & 1 deletion scripts/other/ohc_comparison.R
Original file line number Diff line number Diff line change
Expand Up @@ -153,12 +153,14 @@ obs_data <- get_ohc_data(OHC_FILE, include_unc = T)
default_data <- calc_ohc(NULL, NULL, include_unc= T)
reg_box_data <- calc_ohc(PARAMS, c(0.57, 1.76, 2.38, 2.95, 0.49), include_unc= T)
low_diff_data <- calc_ohc(PARAMS, c(0.57, 1.76, 1.1, 2.95, 0.49), include_unc= T)
ohc_optim_data <- calc_ohc(PARAMS, c(0.65, 1.76, 1.04, 2.33, 0.44), include_unc= T)

default_data$scenario <- "Hector - Default"
reg_box_data$scenario <- "Hector - All Params, Reg Box"
low_diff_data$scenario <- "Hector - All Params, Reg Box, Diff = 1.1"
ohc_optim_data$scenario <- "Hector - All Params, Reg Box, Matilda Diff, Optimize for OHC"

comb_data <- rbind(obs_data, default_data, reg_box_data, low_diff_data)
comb_data <- rbind(obs_data, default_data, reg_box_data, low_diff_data, ohc_optim_data)

ggplot(data = comb_data, aes(x = year, y = value, color = scenario)) +
geom_ribbon(data =
Expand Down

0 comments on commit be355ee

Please sign in to comment.