LUC Reconstruction Analysis in Hector
Pablo Tomberlin
2025-07-24
Hello! Here’s my land-use change (LUC) emissions analysis project using Hector, a simple earth system model. Follow along to see my work.
If you want to run this code yourself, you can clone the repostitory: https://github.com/PabloTomberlin/LUC_Emissions_Analysis
Read in data and load necessary packages
data <- read.csv("selected_GCP.csv")
data_FF <- read.csv("Fossil_Fuels.csv")
observations <- read.csv("GML_observations.csv")
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.2 ✔ tibble 3.2.1
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.0.4
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(Metrics)
library(ggsci)
library(ggrepel)## Warning: package 'ggrepel' was built under R version 4.5.1data <- data %>% filter(!(is.na(year)))Fossil fuel visualization of a lengthening time series (2014, 2019, 2024)
data_FF$GCP <- as.factor(data_FF$GCP)
ggplot() +
geom_line(data = data_FF,
aes(year, value, color = GCP),
linewidth = 0.8) +
facet_wrap("GCP", scales = "fixed") +
scale_color_discrete() +
ggtitle("Fossil Fuel Emissions") +
xlab("Year") +
ylab("Billion tonnes C (GtC/yr)") +
theme_bw() +
theme(text = element_text(size = 12, family = "mono", face = "bold"))
Plot LUC emissions data for 3 key GCP releases (first year, historical extension, and most recent year)
LUC_results <- data %>% filter(GCP %in% c(2007, 2015, 2024))
LUC_results$GCP <- as.factor(LUC_results$GCP)
ggplot() +
geom_line(data = LUC_results,
aes(year, value, group = GCP, color = GCP),
linewidth = 0.8) +
ggtitle("LUC Emissions") +
xlab("Year") +
ylab("Billion tonnes C (GtC/yr)") +
theme_bw() +
theme(text = element_text(size = 12, family = "mono", face = "bold")) +
scale_color_tron()
Basic Hector run
library(hector)
ini_file <- system.file("input/hector_ssp245.ini", package = "hector")
core <- newcore(ini_file)
run(core)## Hector core: Unnamed Hector core
## Start date: 1745
## End date: 2300
## Current date: 2300
## Input file: C:\Users\tomb621\AppData\Local\R\win-library\4.5\hector\input\hector_ssp245.iniout_default <- fetchvars(core, 1850:2100, vars = c(CONCENTRATIONS_CO2(),
GLOBAL_TAS(),
VEG_C(),
SOIL_C(),
LUC_EMISSIONS()))
out_default$GCP <- "Hector Default"LUC plot with Hector Default (SSP 2-45) as separate line
LUC_default <- out_default %>%
filter(variable == "luc_emissions")
ggplot() +
geom_line(data = LUC_results,
aes(year, value, group = GCP),
color = "gray",
linewidth = 0.8) +
geom_label_repel(aes(label = "Hector Unmodified LUC"),
x = 2070,
y = 0.1,
color = "blue") +
facet_wrap("variable", scales = "free") +
geom_line(data = LUC_default,
aes(year, value, group = GCP),
color = "blue",
linewidth = 0.8) +
ggtitle("LUC Emissions") +
xlab("Year") +
ylab("Billion tonnes C (GtC/yr)") +
theme_bw() +
theme(text = element_text(size = 12, family = "mono", face = "bold"))
Create list and for loop to run all carbon budgets individually through Hector and extract projections for climate and carbon cycle variables and conditions
all_gcp_years <- unique(data$GCP)
results <- list()
for(gcp_year in all_gcp_years){
isolated_data <- data %>% filter(GCP == gcp_year)
message("running ", gcp_year)
message("I see ", nrow(isolated_data), " rows of data")
setvar(core, isolated_data$year, LUC_EMISSIONS(), isolated_data$value,
getunits(LUC_EMISSIONS()))
reset(core)
run(core)
results[[as.character(gcp_year)]] <- fetchvars(core, 1850:2100, vars = c(CONCENTRATIONS_CO2(),
GLOBAL_TAS(),
VEG_C(),
SOIL_C(),
LUC_EMISSIONS(),
OCEAN_C()))
}## running 2007## I see 48 rows of data## running 2010## I see 51 rows of data## running 2015## I see 161 rows of data## running 2020## I see 170 rows of data## running 2023## I see 173 rows of data## running 2024## I see 174 rows of dataCombine all data frames including the default model run into one large data frame
results_gcp <- bind_rows(results, .id = "GCP")
results_df <- bind_rows(results_gcp, out_default)
head(results_df)## GCP scenario year variable value units
## 1 2007 Unnamed Hector core 1850 CO2_concentration 283.2175 ppmv CO2
## 2 2007 Unnamed Hector core 1851 CO2_concentration 283.5029 ppmv CO2
## 3 2007 Unnamed Hector core 1852 CO2_concentration 283.7846 ppmv CO2
## 4 2007 Unnamed Hector core 1853 CO2_concentration 284.0911 ppmv CO2
## 5 2007 Unnamed Hector core 1854 CO2_concentration 284.3687 ppmv CO2
## 6 2007 Unnamed Hector core 1855 CO2_concentration 284.6286 ppmv CO2Plot global temperature projections
projections_gcp <- results_df %>%
filter(variable == "global_tas") %>%
filter(GCP %in% c(2007, 2015, 2024))
projections_default <- results_df %>%
filter(variable == "global_tas") %>%
filter(GCP == "Hector Default")
ggplot() +
geom_line(data = projections_gcp,
aes(year, value, color = GCP),
linewidth = 0.8) +
#facet_wrap("variable", scales = "free") +
#geom_line(data = projections_default,
#aes(year, value),
#color = "pink",
#linewidth = 0.8) +
ggtitle("Global temperatures projections") +
xlab("Year") +
ylab("Degrees C") +
theme_bw() +
theme(text = element_text(size = 12, family = "mono", face = "bold")) +
scale_color_tron()
Calculate projection uncertainty due to LUC emissions revisions and store in hector_uncertainty
hector_uncertainty <- results_df %>%
filter(year %in% c(2024, 2100)) %>%
filter(variable %in% c("global_tas", "CO2_concentration")) %>%
group_by(year, variable) %>%
summarise(
mean_value = mean(value),
sd_value = sd(value),
margin_error = 1.96 * (sd_value/sqrt(length(unique(results_df$GCP)))),
spread = max(value) - min(value),
.groups = "drop"
)
print(hector_uncertainty)## # A tibble: 4 × 6
## year variable mean_value sd_value margin_error spread
## <dbl> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 2024 CO2_concentration 433. 7.14 5.29 15.3
## 2 2024 global_tas 1.30 0.0517 0.0383 0.114
## 3 2100 CO2_concentration 558. 13.3 9.85 28.6
## 4 2100 global_tas 2.79 0.0826 0.0612 0.177Relative graph showing deviation each projection has from the average of all temperature projections and a black line to indicate where observations end
results_temp <- results_df %>%
filter(variable == "global_tas")
results_avg <- results_temp %>%
group_by(year) %>%
summarize(mean_value = mean(value))
full_df <- left_join(results_temp, results_avg, by = "year")
diff <- mutate(full_df, subtract = value - mean_value)
tail(diff)## GCP scenario year variable value units
## 1752 Hector Default Unnamed Hector core 2095 global_tas 2.701186 degC
## 1753 Hector Default Unnamed Hector core 2096 global_tas 2.706050 degC
## 1754 Hector Default Unnamed Hector core 2097 global_tas 2.710481 degC
## 1755 Hector Default Unnamed Hector core 2098 global_tas 2.714499 degC
## 1756 Hector Default Unnamed Hector core 2099 global_tas 2.718119 degC
## 1757 Hector Default Unnamed Hector core 2100 global_tas 2.721352 degC
## mean_value subtract
## 1752 2.766486 -0.06529965
## 1753 2.771792 -0.06574168
## 1754 2.776669 -0.06618794
## 1755 2.781137 -0.06663729
## 1756 2.785211 -0.06709258
## 1757 2.788905 -0.06755250ggplot() +
geom_line(data = diff,
aes(year, subtract, color = GCP),
linewidth = 0.8) +
geom_vline(xintercept = 2023, linetype = "dashed", color = "black", linewidth = 0.8) +
geom_hline(yintercept = 0, linetype = "dashed", color = "yellow", linewidth = 0.8) +
xlab("Year") +
ylab("Degrees (C)") +
ggtitle("Difference from average projected temperature") +
theme_bw() +
theme(text = element_text(size = 12, family = "sans", face = "bold"))
Ocean, soil, and vegetation carbon plots to determine the driver of long-lasting LUC effects
soil_veg <- results_df %>%
filter(variable %in% c("soil_c", "veg_c", "ocean_c"))
ggplot() +
geom_line(data = soil_veg,
aes(year, value, color = GCP),
linewidth = 0.8) +
facet_wrap("variable", scales = "free") +
ggtitle("Soil and vegetation carbon projections") +
xlab("Year") +
ylab("Pg C") +
theme_bw() +
theme(text = element_text(size = 12, family = "mono", face = "bold"))
Matilda Analysis (parametric uncertainty) and calculating 95% confidence interval for 100 Matilda runs. This was used to compare how uncertainty of historical LUC emissions and uncertainty of parameters each contributes to overall CO2 projection uncertainty
library(matilda)
ini_file <- system.file("input/hector_ssp245.ini", package = "hector")
core <- newcore(ini_file)
param_sets <- generate_params(core, draws = 100)
m_results <- iterate_model(core,
params = param_sets,
save_vars = c("CO2_concentration"), save_years = 1850:2100)head(m_results)## scenario year variable value units run_number
## 1 Unnamed Hector core 1850 CO2_concentration 285.0622 ppmv CO2 1
## 2 Unnamed Hector core 1851 CO2_concentration 285.3330 ppmv CO2 1
## 3 Unnamed Hector core 1852 CO2_concentration 285.6010 ppmv CO2 1
## 4 Unnamed Hector core 1853 CO2_concentration 285.8958 ppmv CO2 1
## 5 Unnamed Hector core 1854 CO2_concentration 286.1759 ppmv CO2 1
## 6 Unnamed Hector core 1855 CO2_concentration 286.4453 ppmv CO2 1ggplot() +
geom_line(data = m_results,
aes(year, value, group = run_number)) +
xlab("Year") +
ylab("Value (ppmv CO2)") +
ggtitle("Parametric Uncertainty") +
theme_bw() +
theme(text = element_text(size = 12, family = "mono", face = "bold"))
matilda_uncertainty <- m_results %>%
filter(year %in% c(2024, 2100)) %>%
group_by(year, variable) %>%
summarise(
mean_value = mean(value),
sd_value = sd(value),
margin_error = 1.96 * (sd_value/10),
spread = max(value) - min(value),
.groups = "drop"
)
print(matilda_uncertainty)## # A tibble: 2 × 6
## year variable mean_value sd_value margin_error spread
## <dbl> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 2024 CO2_concentration 426. 20.7 4.06 94.9
## 2 2100 CO2_concentration 543. 55.3 10.8 272.Reproducibility
sessionInfo()## R version 4.5.0 (2025-04-11 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26100)
##
## Matrix products: default
## LAPACK version 3.12.1
##
## locale:
## [1] LC_COLLATE=English_United States.utf8
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## time zone: America/New_York
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] matilda_0.1.0 hector_3.2.0 ggrepel_0.9.6 ggsci_3.2.0
## [5] Metrics_0.1.4 lubridate_1.9.4 forcats_1.0.0 stringr_1.5.1
## [9] dplyr_1.1.4 purrr_1.0.4 readr_2.1.5 tidyr_1.3.1
## [13] tibble_3.2.1 ggplot2_3.5.2 tidyverse_2.0.0
##
## loaded via a namespace (and not attached):
## [1] gtable_0.3.6 jsonlite_2.0.0 compiler_4.5.0 Rcpp_1.0.14
## [5] tidyselect_1.2.1 jquerylib_0.1.4 scales_1.4.0 yaml_2.3.10
## [9] fastmap_1.2.0 R6_2.6.1 labeling_0.4.3 generics_0.1.4
## [13] knitr_1.50 bslib_0.9.0 pillar_1.10.2 RColorBrewer_1.1-3
## [17] tzdb_0.5.0 rlang_1.1.6 utf8_1.2.5 cachem_1.1.0
## [21] stringi_1.8.7 xfun_0.52 sass_0.4.10 timechange_0.3.0
## [25] cli_3.6.5 withr_3.0.2 magrittr_2.0.3 digest_0.6.37
## [29] grid_4.5.0 rstudioapi_0.17.1 hms_1.1.3 lifecycle_1.0.4
## [33] vctrs_0.6.5 evaluate_1.0.3 glue_1.8.0 farver_2.1.2
## [37] rmarkdown_2.29 tools_4.5.0 pkgconfig_2.0.3 htmltools_0.5.8.1