Does heat flux solve the S and diff problem?

Purpose

Now that we have determined that the land an ocean tempearture does not resolve the issue of S and diff we are going to try using heat flux as comparison data.

Set Up and Define Functions

library(hector)
devtools::load_all("/Users/dorh012/Documents/2019/hectorcal")

Calibrate Hector to CESM1-CAM5 data.

# To make things easy we are going to start out just looking at the rcp 85 + historical results. 
cmip_individual %>% 
    filter(model == 'CESM1-CAM5' & experiment %in% c('historical', 'rcp85')) %>%  
    select(year, model, value) %>% 
    mutate(variable = GLOBAL_TEMP(), experiment = 'rcp85') -> 
    cesm_comp_data
core <- newcore(system.file('input/hector_rcp85.ini', package = 'hector'))
fn   <- minimize_global(core = core, comparison_data = cesm_comp_data)
params <- c(3.5, 1.5, 0.6)
names(params) <- c(ECS(), AERO_SCALE(), VOLCANIC_SCALE())
print(params)

     S  alpha volscl 
   3.5    1.5    0.6 

cesm_temp_fits <- optim(par = params, fn = fn)

What were the best fit params?

cesm_temp_fits$par

       S    alpha   volscl 
3.649249 1.221565 1.033603 

Emulate CESM1-CAM5 Heat Flux Data

Use the best fits to run Hector again and pull out the global temp and the heat flux.

parameterize_core(params = cesm_temp_fits$par, core)

Hector core:    unnamed hector core
Start date: 1745
End date:   2300
Current date:   1745
Input file: /Users/dorh012/Library/R/3.5/library/hector/input/hector_rcp85.ini

reset(core)

Hector core:    unnamed hector core
Start date: 1745
End date:   2300
Current date:   1745
Input file: /Users/dorh012/Library/R/3.5/library/hector/input/hector_rcp85.ini

run(core)

Hector core:    unnamed hector core
Start date: 1745
End date:   2300
Current date:   2300
Input file: /Users/dorh012/Library/R/3.5/library/hector/input/hector_rcp85.ini

fetchvars(core, 1850:2100, c(GLOBAL_TEMP(), HEAT_FLUX())) %>%
    select(year, value, variable) %>%
    mutate(experiment = 'rcp85', model = 'emulated CESM') ->
    emulatedCESM_comp_data

Fit S, alpha, and volscl to CESM1-CAM5 at fixed diff

diff_values <- seq(from = 0.5, to = 15, by = 1)
lapply(diff_values, function(x){
    # Make a new core
    core <- newcore(system.file('input/hector_rcp85.ini', package = 'hector'))
    # Reset the diff value
    dif        <- x
    names(dif) <- DIFFUSIVITY()
    parameterize_core(core = core, params = dif)
    reset(core)
    # The function to minimize
    fn <- minimize_temp_flux(core = core, comparison_data = emulatedCESM_comp_data)
    # Make the inital guess for the parameter
    best_guess        <- c(3.649204, 1.221605, 1.033848)
    names(best_guess) <- c( 'S', 'alpha', 'volscl')
    # Minmize the paramters
    fit_heatFlux_temp <- optim(par = best_guess, fn = fn,  control = list('maxit' = 900))
    # If the optimization was sucessful then use it as hector input.
    if(fit_heatFlux_temp$convergence == 0 ){
        parameterize_core(params = fit_heatFlux_temp$par, core)
        reset(core)
        run(core, runtodate = 2100)
        fetchvars(core, dates = 1750:2100, vars = c(hector::GLOBAL_TEMP(), hector::HEAT_FLUX())) %>%
            mutate(S = fit_heatFlux_temp$par[['S']],
                   diff = x,
                   alpha = fit_heatFlux_temp$par[['alpha']],
                   volscl = fit_heatFlux_temp$par[['volscl']],
                   SE = fit_heatFlux_temp$value)
    }}) %>%
    bind_rows ->
    heatFlux_temp_fits

How do the fits look?

ggplot(data = heatFlux_temp_fits) +
    geom_point(data = cesm_comp_data, aes(year, value, color = 'actual CESM1')) +
    geom_point(data = emulatedCESM_comp_data, aes(year, value, color = 'emulated CESM1')) +
    geom_line(aes(year, value, group = diff)) +
    facet_wrap('variable') +
    coord_cartesian(xlim = c(1850, 2100)) +
    labs(title = 'Hector fits at fixed diff vs comparison data')

What is the relationship between S and diff?

heatFlux_temp_fits %>%
    select(S, diff, 'minimized squared error' = SE) %>%
    distinct %>%
    ggplot(aes(diff, S, color = `minimized squared error`)) +
    geom_point(size = 3) +
    coord_cartesian(ylim = c(0, 6)) +
    theme(legend.position = 'bottom') +
    labs(title = 'Best fit S vs diff')

The darker the blue the better the fit is. The dark blue dots are over the range of diff values that we expect to be reasonable.

What happens if we fit for alpha and S?

# Inital parameter guesses
inital_guess        <- c(3.5, 1.6,  0.6, 2.3)
names(inital_guess) <- c(ECS(), AERO_SCALE(), VOLCANIC_SCALE(), DIFFUSIVITY())
# Make the function to minimize and optmize the fit.
fn_heatF_tmep <- minimize_temp_flux(core = core, comparison_data = emulatedCESM_comp_data)
fit_S_diff    <- optim(par = inital_guess, fn = fn_heatF_tmep)

What do the parameter values look like?

fit_S_diff$par

       S    alpha   volscl     diff 
3.649284 1.221651 1.033753 2.299885 

Other than the alpha value, the other parameter values are pretty realistic.

What does the fit look like?

# Run hector with parameter values
parameterize_core(fit_S_diff$par, core)

Hector core:    unnamed hector core
Start date: 1745
End date:   2300
Current date:   1745
Input file: /Users/dorh012/Library/R/3.5/library/hector/input/hector_rcp85.ini

run(core)

Hector core:    unnamed hector core
Start date: 1745
End date:   2300
Current date:   2300
Input file: /Users/dorh012/Library/R/3.5/library/hector/input/hector_rcp85.ini

hector_output <- fetchvars(core = core, dates = 1850:2100, vars = c(GLOBAL_TEMP()))
ggplot(data = hector_output, aes(year, value, color = 'hector output')) +
    geom_point(data = cesm_comp_data, aes(year, value, color = 'CESM1-CAM5')) +
    geom_line(size = 1)

Looks pretty good to me!

Conclusions

Using heat flux as comparison data resolved the problem of S and diff! whoot whoot!

Next Steps

1. Get heat flux scale and center results? (do we want it from the PCA or I could just get it from the ensemble)
2. Test with the actual CESM1-CAM5 data
3. Download/process the CMIP5 results for the heat flux data