library(knitr)
library(here)
library(dplyr)
library(ggplot2)
library(hector)

opts_chunk$set(
  collapse = TRUE,
  warning = FALSE,
  message = FALSE,
  echo = FALSE,
  comment = "#>"
)

pkgload::load_all(here())

fig <- function(...) {
  filepath <- here("analysis", "figures", ...)
  knitr::include_graphics(file_in(filepath))
}

Introduction

Permafrost---soil that continuously remains below 0°C for at least two consecutive years---underlies an area of 22 (± 3) million km~2~, roughly 17% of the Earth's exposed land surface [@gruber_2012_derivation], and contains an estimated 1300 (1100-1500) Pg of organic carbon [@hugelius_2014_estimated]. Recent increases in global air temperature [@ipcc_ar5_wg1], which are amplified at high latitudes [@pithan_2014_arctic], have resulted in widespread permafrost thaw [@romanovsky_2010_permafrost], and simulations from variety of climate and land surface models across a wide range of scenarios suggest that this trend will continue into the future [@koven_2013_analysis; @chadburn_2017_observation]. Permafrost thaw can dramatically alter surface terrain and hydrology [@jones_2013_quantifying; @godin_2014_effects; @necsoiu_2013_multitemporal], with adverse consequences for human infrastructure in permafrost regions [@anisimov_1997_permafrost; @nelson_2002_climate; @larsen_2008_estimating]. Moreover, as permafrost thaws, this carbon becomes available to microbes for decomposition, resulting in the production of CO~2~ and CH~4~ [@brown_2008_report; @romanovsky_2010_thermal; @bond-lamberty_2016_temperature] that could lead to further warming [@koven_2011_permafrost; @schuur_2015_climate]. Accounting for this permafrost carbon-climate feedback generally increases projections of greenhouse gas concentrations and global temperatures (REFS) and increases estimates of the economic impact of climate change [@hope_2015_economic; @yumashev_2019_climate; @chen_2019_economic].

The higher carbon emissions associated with warming-driven permafrost thaw may be offset by increases in primary productivity. Some studies based on meteorological tower measurements of carbon flux and optical remote sensing suggest that high-latitude ecosystems (mainly tundra and boreal forest) remain carbon-neutral or are even a small net carbon sink [@mcguire_2012_assessment; @welp_2016_increasing]. However, other studies, based on different sets of flux measurements and airborne gas sampling, suggest that the high-latitude regions are a net carbon source [@belshe_2013_tundra; @commane_2017_carbon; @natali_2019_large]. The uncertainty in the net high-latitude carbon flux may be driven by the heterogeneity of high-latitude landscapes in terms of vegetation cover, soil properties, topography, and many other factors known to affect both the rate of permafrost thaw and the subsequent carbon flux [@turetsky_2002_boreal; @wickland_2006_effects; @lund_2010_variability; @james_2013_multidecadal; @johnson_2013_permafrost; @grant_2019_modeling1; @grant_2019_modeling2]. This uncertainty is also present in recent permafrost modeling studies [@burke_2017_quantifying; @qian_2010_enhanced; @ito_2016_impacts; @harp_2016_effect].

Land surface models are expensive to run, making it challenging to use them for uncertainty quantification and exploration of alternative policy scenarios. Simple climate models are an alternative. (More on simple climate models). (More on permafrost in simple climate models).

Hector [@hartin_2015_simple; @vega-westhoff_2019_impacts]. In this study, we implement a simple representation of permafrost thaw and evaluate its consequences for global climate in Hector.

Methods

Model description: Hector

Hector [@hartin_2015_simple,@hartin_2016_ocean,]. Simple climate model.

(TODO: More details on terrestrial C cycle in Hector). The default heterotrophic respiration ($R$) scheme for a pool $p$ (detritus or soil) in Hector:

$$ R_{p} = C_p \times f_p \times Q_{10} ^ \frac{T}{10} $$

Hector permafrost sub-model

The general principle is that permafrost constitutes an additional reserve of soil carbon that, because it is frozen, is inaccessible to microbes. As permafrost thaws and some fraction of this carbon becomes accessible to microbes, it is transferred from the frozen permafrost pool to the standard soil C pool, where it is decomposed by Hector's standard decomposition routine. Let $C_{pf}[t]$ be the permafrost C pool and $C_{s}[t]$ be the soil C pool at time $t$ (both in units of Pg C), and let $\Delta C_{pf}[t]$ be the change in the permafrost C pool at time $t$. The C consequences of permafrost thaw can therefore be represented as:

[ C_{pf}[t] = C_{pf}[t-1] - \Delta C_{pf}[t] ]

[ C_{s}[t] = C_{s}[t] + \Delta C_{pf}[t] ]

Let $C[0]{pf}$ be the initial size of the permafrost C pool (Pg C) and $\Phi[t]$ be the fraction of permafrost remaining (in arbitrary area or volume units) at time $t$. Assuming a uniform permafrost C density, $\Delta C{pf}[t]$ can be expressed as:

[ \Delta C_{pf}[t] = \Delta \Phi[t] C_{pf}[t-1] ]

To a first approximation, $\Phi[t]$ is a function of air temperature ($T[t]_{air}$). Kessler (2015) assume this relationship is linear. However, because the permafrost area fraction is, by definition, bounded by zero (and 1 [^1]), and because deeper permafrost thaws more slowly than shallow permafrost, we use the lognormal cumulative distribution function ($NCDF(\log(x) | \mu, \sigma)$) instead:

[^1]: Technically, permafrost area could increase in the case of cooling temperatures, and therefore the area fraction could be greater than 1. However, because even the most aggressive climate action scenarios show temperatures that stabilize above year 2000, we assume that permafrost area will never grow more than the starting value.

[ \Phi[t] = 1 - NCDF(\log(x) | \mu, \sigma) ]

where $\mu$ and $\sigma$ are model parameters whose values are described in the Configuration section below. The change in frozen fraction at a given time step, $\Delta \Phi[t]$, is given by:

[ \Delta \Phi[t] = \Phi[t] - \Phi[t-1] ]

Methane emissions

In the default configuration of Hector, natural CH~4~ emissions are prescribed at a constant value (300 Tg year^-1^) [@hartin_2015_simple]. In the new configuration, CH~4~ emissions are a fixed fraction of heterotrophic respiration.

Total heterotrophic respiration ($RH[i,t]$) for biome $i$ at time $t$ is the sum of heterotrophic respiration in detritus ($RH_d$) and soil ($RH_s$).

$$ RH[i,t] = RH_s[i,t] + RH_d[i,t] $$

Detritus and soil heterotrophic respiration are proportional to the sizes of their respective carbon pools ($C_d$ and $C_s$, respectively, both in Pg C), with a rate that increases exponentially with temperature according to a biome-specific temperature sensitivity parameter ($Q_{10}$). Detritus respiration increases with biome-specific air temperature change since pre-industrial ($T[i,t]$), while soil respiraiton increases with the 200 year running mean of air temperature ($T_{200}[i,t]$) because soils warm slowly relative to the atmosphere.

[ T_{air}[i,t] = \delta \, T_{air}[t] ]

[ RH_d[i,t] = \frac{1}{4} C_d Q_{10}[i] ^ {{T}[i,t] / 10} ]

[ RH_s[i,t] = \frac{1}{50} C_s Q_{10}[i] ^ {T_{200}[i,t] / 10} ]

In the original formulation, natural CH~4~ emissions were a prescribed constant value. In the revised version here, natural CH~4~ emissions are a fraction of heterotrophic respiration ($RH_{CH_4}[i,t]$) based on the biome-specific CH~4~ fraction parameter ($f_{CH_4}[i]$):

[ RH_{CH_4}[i,t] = f_{CH_4}[i] RH[i,t] ]

Configuration

Initial permafrost C is set to 1035 Pg C based on the estimate for near-surface (<3 m depth) permafrost by @hugelius_2014_estimated.

# Permafrost exponential model
pfa <- 2.371
pfb <- -0.676
pfc <- -3.685

# Kessler linear model
kessler_m <- 0.172
kessler_b <- 0.8 * kessler_m + 1

For globally-averaged permafrost, exponential parameters $a = r pfa$, $b = r pfb$, and $c = r pfc$ to most closely reproduce the rate of r kessler_m K^-1^ reported by Kessler (2015) over the range of 0.8 [^2] to 4 K above the pre-industrial baseline.

[^2]: Kessler (2015) report this as temperature change from year 2000. 0.8 K is the warming since pre-industrial as estimated by the default Hector configuration.

sigmoid_cap <- paste(
  "Fraction of permafrost thaw as a function of change in",
  "global annual mean air temperature since pre-industrial (1750)."
)
plot(0, 0, type = "n", xlim = c(0, 8), ylim = c(0, 1),
     xlab = expression(Delta * T[air] ~ (K)),
     ylab = "Frac. permafrost remaining")
curve(1 - (1 + pfa * exp(pfb * x)) ^ pfc, 0, 10, col = 1, add = TRUE)
abline(a = kessler_b, b = -kessler_m, lty = "dashed")
abline(h = 1, lty = "dotted")
legend("topright", c("Hector", "Kessler (2015)"), lty = c("solid","dashed"),
       bg = "white")

Results

Effect of permafrost C

permafrost_c_cap <- paste(
  "Effect of permafrost C emissions on scenarios."
)

run_rcp <- function(rcp) {
  inidir <- system.file("input", package = "hector")
  inifile <- normalizePath(file.path(inidir, paste0("hector_rcp", rcp, ".ini")))
  stopifnot(file.exists(inifile))
  hc <- newcore(inifile)
  dates <- seq(1750, 2100)
  run(hc, max(dates))
  outvars <- c(ATMOSPHERIC_CO2(), GLOBAL_TEMP(), F_FROZEN(), SOIL_C())
  no_pf <- fetchvars(hc, dates, outvars, scenario = "Original")
  reset(hc)
  setvar(hc, 0, PERMAFROST_C(), 1035, "PgC")
  run(hc, max(dates))
  yes_pf <- fetchvars(hc, dates, outvars, scenario = "Permafrost")
  return(tibble::as_tibble(dplyr::bind_rows(no_pf, yes_pf)))
}

rcps <- c("26", "45", "60", "85")
names(rcps) <- paste0("RCP", rcps)
results <- purrr::map_dfr(rcps, run_rcp, .id = "RCP")

results %>%
  mutate(variable = recode_factor(
    variable,
    Tgav = "Global mean temperature anomaly (K)",
    Ca = "Atmospheric CO2 (ppm)",
    soil_c = "Active soil C pool (PgC)",
    permafrost_c = "Permafrost C pool (PgC)",
    f_frozen = "Frozen permafrost fraction"
  )) %>%
  ggplot() +
  aes(x = year, y = value, linetype = scenario, color = RCP) +
  geom_line() +
  facet_wrap(vars(variable), scales = "free_y") +
  theme_bw() +
  theme(axis.title = element_blank())

Discussion

Rate of permafrost C release also depends on soil moisture conditions -- drier soils release C much faster ("carbon bomb") than wetter soils ("carbon fizz") [@elberling_2013_long]. Moisture will also affect the balance of aerobic (CO2 release) vs. anaerobic (CH4) C release [@turetsky_2002_boreal], to the extent that unclear if anaerobic (wet) areas are C sources or sinks [@wickland_2006_effects]. Effects of permafrost thaw on soil moisture are a complex hydrological problem -- drainage very sensitive to local (micro-)topography [@wickland_2006_effects]. So will vegetation cover [@wickland_2006_effects].

Temperature amplification of permafrost carbon feedback (by 2100) 0.02 to 0.36 °C [@burke_2013_estimating; @schneider_2012_estimating; @schneider_2015_observation], or 0.1 to 0.8 °C in [@macdougall_2012_significant; @macdougall_2013_if], or 10-40% of peak temperature change [@crichton_2016_permafrost], or 0.2 to 12% [@burke_2017_quantifying].

Permafrost carbon has greater impact at lower emissions scenarios [@burke_2017_quantifying; @macdougall_2012_significant; @macdougall_2013_if] .

Conclusion

Acknowledgments

Funded by EPA grant XXX. Cyberinfrastructure support from Pacific Northwest National Laboratory (PNNL).

Miscellaneous notes

Permafrost emissions scenarios

Digitized scenarios from [@schaefer_2011_amount]. SiBCASA model predictions of CO2 emissions (permafrost respiration; $R_{pc}$; note -- no methane!) through 2300. These results were digitized using WebPlotDigitizer (https://apps.automeris.io/wpd/), and interpolated to annual resolution (using R stats::spline function).

Digitized scenarios from [@hope_2015_economic]. CO2 and CH4 emissions from SiBCASA model.

[@schuur_2009_effect] -- Estimate 0.8 - 1.1 Pg C yr-1.

Back-of-the-envelope estimates from [@zimov_2006_climate_change]: - 500 Gt C in loess that could be completely emitted by 2100 (plus other C sources). - 10-40 g C m-3 day-1 off the bat, slowing down to equilibrium (?) rate of 0.5-5 g C m-3 day-1 for several years. Combine with data on permafrost spatial extent, density, etc. to generate estimates (but can back-calculate from 500 Gt C above?)

Parameter calibration

Some of these are time series, while others are individual estimates at particular points in time. To give them equal weight in the likelihood, we down-weight the time series likelihoods by the number of time steps.

We derived a distribution for the Arctic warming factor from TODO.

TODO: Table and multi-panel figure of input datasets.

For the $\alpha$ and $\phi$ parameters in case 3, we looked at the literature on permafrost methane emissions [e.g., @wickland_2006_effects].

Other notes

Frozen carbon residence time (FCRt) from [@burke_2017_quantifying]:

$$ FCRt = FCRt0 * exp(-\Delta T / \Gamma) (for \Delta T > 0.2 °C) $$

Other Hector parameters to consider.

| Variable | INI name | Description | Value | |----------|----------|---------------------------------------------|-------| | $f_{nv}$ | f_nppv | Fraction of NPP C transferred to vegetation | 0.35 | | $f_{nd}$ | f_nppd | Fraction of NPP C transferred to detritus | 0.60 | | $f_{nd}$ | | Fraction of NPP C transferred to soil | 0.05 | | $f_{lv}$ | f_lucv | Fraction of LUC change flux from vegetation | 0.10 | | $f_{ld}$ | f_lucd | Fraction of LUC change flux from detritus | 0.01 | | $f_{ls}$ | | Fraction of LUC change flux from soil | 0.89 | | $f_{ds}$ | | Fraction of detritus C that goes to soil | 0.60 | | $f_{rd}$ | | Fraction of respiration C to detritus | 0.25 | | $f_{rs}$ | | Fraction of respriation C to soil | 0.02 |

According to [@hartin_2015_simple], these were selected to be "generally consistent with previous simple earth system models [e.g., @meinshausen_2011_emulating]".

pagebreak

References

pagebreak

Colophon

This report was generated on r Sys.time() using the following computational environment and dependencies:

# which R packages and versions?
devtools::session_info()

The current Git commit details are:

# what commit is this file at? You may need to change the path value
# if your Rmd is not in analysis/paper/
git2r::repository("../..")


ashiklom/hector_permafrost_emit documentation built on March 26, 2020, 12:15 a.m.