Identical outputs for climate change scenarios using CMIP6 data
Francisco J. Muñoz Nolasco
I been trying to use NicheMapR to model the thermoregulation and water balance of the Mexican salamander Isthmura belli under current and future climate scenarios. The current scenario model runs successfully using NCEP data via micro_ncep.
For the future scenarios, my initial approach was to download CMIP6 data for various SSPs and time periods using geodata::cmip6_tile. I then calculated the mean temperature delta between the future CMIP6 data and a historical baseline. I attempted to apply this delta to the temperature values within the micro object generated by micro_ncep before running the ectotherm function.
However, all my future scenario model runs (for SSP1-2.6, SSP3-7.0, and SSP5-8.5 across three time periods: 2041-2060, 2061-2080, and 2081-2100) are returning identical outputs. Upon debugging, it appears that the micro object from the current scenario is not being correctly modified for each iteration of the loop, leading to the same input for every future run.
The issue is that CMIP6 data is static, representing long-term means. It cannot be directly used to generate a daily-time step micro object, which is required by ectotherm. As far as I know, the NicheMapR documentation and common practices suggest using the micro_clim function to generate a full microclimate dataset based on a time series of daily weather variables, but I'm not sure how to feed the CMIP6 climate change anomalies (the delta values) into a daily time series like the NCEP data.
What is the correct and most robust method to project the CMIP6 anomalies onto a NCEP-like daily time series for use in NicheMapR? Specifically, how should I properly update the micro object within a loop to reflect the different future climate scenarios, including changes in temperature, precipitation, and solar radiation, while maintaining the daily variability? Any guidance on the recommended workflow for creating future microclimate inputs would be greatly appreciated. I attach my code. Thank you in advance.
# Load packages
library(NicheMapR)
library(geodata)
library(terra)
# Current scenario
loc <- c(-103.9456, 19.4453)
dstart <- "01/01/2023"
dfinish <- "31/12/2023"
altt <- 2440
refl <- 0.10
slr <- 0.96
runshade <- 1
runmoist <- 1
dep <- c(0, 2.5, 5, 10, 15, 20, 30, 50, 100, 200)
soiltype <- 4
params <- list(
ID = "salamander_default",
shape = 4,
mass = 0.029337,
pct_wet = 0.85,
A_max = 0.9,
A_min = 0.9,
Z_b = 0.02,
sub_resist = 1,
T_B_min = 10.2,
T_F_min = 10.2,
T_F_max = 25.3,
T_pref = 19.34,
CT_max = 30,
CT_min = 4.78,
mindepth = 0.01,
maxdepth = 0.80,
shade_seek = 1,
burrow = 1,
climb = 0,
nocturn = 1,
crepus = 1,
diurn = 0
)
# Run current model
micro_current <- micro_ncep(loc = loc, dstart = dstart, dfinish = dfinish)
salamander_current <- ectotherm(micro = micro_current, params = params)
salamander_current_data <- as.data.frame(salamander_current$environ)
# Future scenarios setup
ssps <- c("126", "370", "585")
time_periods <- c("2041-2060", "2061-2080", "2081-2100")
future_results <- data.frame(
scenario = character(),
time_period = character(),
mean_body_temp = numeric(),
total_active_hours = numeric(),
average_daily_water_loss = numeric()
)
# Future scenarios loop
for (ssp in ssps) {
for (time_period in time_periods) {
message(paste0("Running SSP ", ssp, " period ", time_period, " ..."))
# Download future CMIP6 data
cmip_future <- cmip6_tile(
lon = loc[1], lat = loc[2],
model = "CNRM-CM6-1", ssp = ssp,
time = time_period, var = c("tmax", "tmin", "prec", "srad", "wind"), path = tempdir()
)
# Download baseline CMIP6 data (same time period as NCEP)
cmip_baseline <- cmip6_tile(
lon = loc[1], lat = loc[2],
model = "CNRM-CM6-1", ssp = ssp,
time = "2021-2040", var = c("tmax", "tmin", "prec", "srad", "wind"), path = tempdir()
)
# Calculate daily weather deltas from CMIP6 data
temp_delta_vec <- global(cmip_future$tmax, mean, na.rm = TRUE) - global(cmip_baseline$tmax, mean, na.rm = TRUE)
prec_delta_vec <- global(cmip_future$prec, mean, na.rm = TRUE) - global(cmip_baseline$prec, mean, na.rm = TRUE)
srad_delta_vec <- global(cmip_future$srad, mean, na.rm = TRUE) / global(cmip_baseline$srad, mean, na.rm = TRUE)
# Apply deltas to NCEP data
micro_future <- micro_current
micro_future$metout[, 'TAREF'] <- micro_current$metout[, 'TAREF'] + as.vector(temp_delta_vec)
micro_future$metout[, 'PCT_PREC'] <- micro_current$metout[, 'PCT_PREC'] + as.vector(prec_delta_vec)
micro_future$metout[, 'SRAD'] <- micro_current$metout[, 'SRAD'] * as.vector(srad_delta_vec)
# Run ectotherm model for future scenario
salamander_future <- ectotherm(micro = micro_future, params = params)
# Extract and store results
future_data <- as.data.frame(salamander_future$environ)
mean_temp <- mean(future_data$TC, na.rm = TRUE)
active_hours <- sum(future_data$ACT > 0, na.rm = TRUE)
mass_bal_future <- as.data.frame(salamander_future$masbal)
mass_bal_future$H2O_total <- mass_bal_future$H2OCut_g + mass_bal_future$H2OResp_g
daily_water_loss_future <- aggregate(H2O_total ~ as.Date(micro_future$dates), data = mass_bal_future, sum)
avg_water_loss <- mean(daily_water_loss_future$H2O_total, na.rm = TRUE)
future_results <- rbind(future_results, data.frame(
scenario = paste0("SSP", ssp),
time_period = time_period,
mean_body_temp = mean_temp,
total_active_hours = active_hours,
average_daily_water_loss = avg_water_loss
))
}
}
# Check results
print(future_results)
NicheMapR
TAIRhr_orig <- TAIRhr
yearstodo <- seq(ystart, yfinish)
nyears <- yfinish - ystart + 1
if(yfinish > 2015){
ystart_terra <- 2015 - nyears + 1
yfinish_terra <- 2015
message(paste0("terraclimate climate change scenarios are for 1985 to 2015 - using ", ystart, "-", yfinish), '\n')
}else{
ystart_terra <- ystart
yfinish_terra <- yfinish
}
leapyears <- seq(1900, 2060, 4)
message("generate climate change scenario", '\n')
# diff spline function
getdiff <- function(diffs){
diff1 <- (unlist(diffs[1]) + unlist(diffs[length(diffs)])) / 2
# generate list of days
leapcount <- 0
for(ys in 1:nyears){
if(ys == 1){
if(nyears == 1){
if(yearstodo[ys] %in% leapyears){
leapcount <- leapcount + 1
day<-c(1, 15.5, 45.5, 75.5, 106, 136.5, 167, 197.5, 228.5, 259, 289.5, 320, 350.5, 366)
}else{
day<-c(1, 15.5, 45, 74.5, 105, 135.5, 166, 196.5, 227.5, 258, 288.5, 319, 349.5, 365)
}
}else{
if(yearstodo[ys] %in% leapyears){
leapcount <- leapcount + 1
day<-c(1, 15.5, 45.5, 75.5, 106, 136.5, 167, 197.5, 228.5, 259, 289.5, 320, 350.5)
}else{
day<-c(1, 15.5, 45, 74.5, 105, 135.5, 166, 196.5, 227.5, 258, 288.5, 319, 349.5)
}
}
}else{
if(ys == nyears){
if(yearstodo[ys] %in% leapyears){
leapcount <- leapcount + 1
day<-c(15.5, 45.5, 75.5, 106, 136.5, 167, 197.5, 228.5, 259, 289.5, 320, 350.5, 366)
}else{
day<-c(15.5, 45, 74.5, 105, 135.5, 166, 196.5, 227.5, 258, 288.5, 319, 349.5, 365)
}
}else{
if(yearstodo[ys] %in% leapyears){
leapcount <- leapcount + 1
day<-c(15.5, 45.5, 75.5, 106, 136.5, 167, 197.5, 228.5, 259, 289.5, 320, 350.5)
}else{
day<-c(15.5, 45, 74.5, 105, 135.5, 166, 196.5, 227.5, 258, 288.5, 319, 349.5)
}
}
}
if(ys == 1){
days2 <- day
days <- day
}else{
if(yearstodo[ys] %in% leapyears){
days2 <- c(days2, (day + 366 * (ys - 1)) + leapcount)
}else{
days2 <- c(days2, (day + 365 * (ys - 1)) + leapcount)
}
days <- c(days, day)
}
}
diffs3 <- c(diff1, diffs, diff1)
days_diffs <- data.frame(matrix(NA, nrow = nyears * 12 + 2, ncol = 3))
days_diffs[, 1] <- days
days_diffs[, 3] <- days2
days_diffs[, 2] <- diffs3
colnames(days_diffs) <- c("days", "diffs", "new_day")
# interpolate monthly differences
f <- approxfun(x = days_diffs$new_day, y = days_diffs$diffs)
xx <- seq(1, max(days2), 1)
sp_diff <- f(xx)
return(sp_diff)
}
terra <- as.data.frame(get_terra(x = loc, ystart = ystart_terra, yfinish = yfinish_terra, source = terra_source))
if(is.infinite(max(terra))){
message("no TerraClimate data available", '\n')
stop()
}
if(scenario == 4){
terra_cc <- as.data.frame(get_terra(x = loc, ystart = ystart_terra, yfinish = yfinish_terra, scenario = 4, source = terra_source))
}
if(scenario == 2){
terra_cc <- as.data.frame(get_terra(x = loc, ystart = ystart_terra, yfinish = yfinish_terra, scenario = 2, source = terra_source))
}
tmin_diffs <- terra_cc$TMINN - terra$TMINN
tmax_diffs <- terra_cc$TMAXX - terra$TMAXX
precip_diffs <- terra_cc$RAINFALL / terra$RAINFALL
srad_diffs <- terra_cc$SRAD / terra$SRAD
vpd_diffs <- terra_cc$VPD / terra$VPD
srad_diffs[is.na(srad_diffs)] <- 1
vpd_diffs[is.na(vpd_diffs)] <- 1
precip_diffs[is.na(precip_diffs)] <- 1
srad_diffs[is.infinite(srad_diffs)] <- 1
vpd_diffs[is.infinite(vpd_diffs)] <- 1
precip_diffs[is.infinite(precip_diffs)] <- 1
TMINN_diff <- getdiff(tmin_diffs)
TMAXX_diff <- getdiff(tmax_diffs)
VPD_diff <- getdiff(vpd_diffs)
SOLAR_diff <- getdiff(srad_diffs)
RAIN_diff <- getdiff(precip_diffs)
TMINN_diff <- rep(TMINN_diff, each=24)[1:length(TAIRhr)]
TMAXX_diff <- rep(TMAXX_diff, each=24)[1:length(TAIRhr)]
VPD_diff <- rep(VPD_diff, each=24)[1:length(TAIRhr)]
SOLAR_diff <- rep(SOLAR_diff, each=24)[1:length(TAIRhr)]
# code to apply changes in min and max air temperature to hourly air temperature data
# find times of maxima and minima
mins <- aggregate(TAIRhr, by = list(format(tt, "%Y-%m-%d")), FUN = min)$x
maxs <- aggregate(TAIRhr, by = list(format(tt, "%Y-%m-%d")), FUN = max)$x
mins <- rep(mins, each=24)
maxs <- rep(maxs, each=24)
mins2 <- TAIRhr / mins[1:length(TAIRhr)]
maxs2 <- TAIRhr / maxs
mins2[mins2 != 1] <- 0
maxs2[maxs2 != 1] <- 0
# construct weightings so that changes in min and max air temp can be applied
minweight <- rep(NA, length(mins2))
maxweight <- minweight
minweight[mins2 == 1] <- 1
minweight[maxs2 == 1] <- 0
maxweight[mins2 == 1] <- 0
maxweight[maxs2 == 1] <- 1
minweight <- zoo::na.approx(minweight, na.rm = FALSE)
maxweight <- zoo::na.approx(maxweight, na.rm = FALSE)
minweight[is.na(minweight)] <- 0.5
maxweight[is.na(maxweight)] <- 0.5
minweight[minweight > 1] <- 1
minweight[minweight < 0] <- 0
maxweight[maxweight > 1] <- 1
maxweight[maxweight < 0] <- 0
# construct final weighted correction factor and apply
diff <- TMINN_diff * minweight + TMAXX_diff * maxweight
TAIRhr <- TAIRhr + diff
}
SOLRhr <- hourlyradwind$swrad / 0.0036
if(scenario != 0){
SOLRhr <- SOLRhr * SOLAR_diff
}
SOLRhr[SOLRhr < 0] <- 0
SOLRhr <- zoo::na.approx(SOLRhr)
IRDhr <- hourlydata$downlong / 0.0036
if(IR != 2){
IRDhr <- IRDhr * 0 + -1 # make negative so computed internally in the microclimate model
}
CLDhr <- hourlydata$cloudcover
if(scenario != 0){
CLDhr <- CLDhr * (1 + 1 - SOLAR_diff)
}
CLDhr[CLDhr < 0] <- 0
CLDhr[CLDhr > 100] <- 100
RHhr <- suppressWarnings(humidityconvert(h = hourlydata$humidity, intype = 'specific', p = hourlydata$pressure, tc = TAIRhr)$relative)
RHhr[RHhr > 100] <- 100
RHhr[RHhr < 0] <- 0
if(scenario != 0){
VPDhr <- WETAIR(db = TAIRhr_orig, rh = 100)$e - WETAIR(db = TAIRhr_orig, rh = RHhr)$e
VPDhr <- VPDhr - mean(VPD_diff) # note using mean here because otherwise can lead to unrealistically low RH which then stronly affects TSKY
e <- WETAIR(db = TAIRhr, rh = 100)$e - VPDhr
es <- WETAIR(db = TAIRhr, rh = 100)$esat
RHhr <- (e / es) * 100
RHhr[RHhr > 100] <- 100
RHhr[RHhr < 0] <- 0
}
WNhr <- hourlyradwind$windspeed
WNhr[is.na(WNhr)] <- 0.1
if(rainhourly == 0){
RAINhr <- rep(0, 24 * ndays)
}else{
RAINhr <- rainhour
}
if(length(rainhourly) == length(tme)){ # an hourly rainfall vector has been provided
aggvec <- rep(1:(length(tme)/24), 24)
aggvec <- aggvec[order(aggvec)]
RAINFALL <- aggregate(rainhourly, by = aggvec, FUN = 'sum') # aggregate to daily totals
}else{
RAINFALL <- dailyprecip
}
PRESShr <- hourlydata$pressure
if(scenario != 0){
RAINFALL <- RAINFALL * RAIN_diff
}
RAINFALL[RAINFALL < 0.1] <- 0