In special-uor/codos: CO2 Correction Tools

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-",
  out.width = "100%",
  dpi = 500,
  dev = "cairo_pdf"
)

CO_dos: CO₂ Correction Tools

r if(!knitr::is_latex_output()) badger::badge_devel("special-uor/codos", "yellow") r if(!knitr::is_latex_output()) badger::badge_cran_release("codos", "red") r if(!knitr::is_latex_output()) badger::badge_github_actions("special-uor/codos") r if(!knitr::is_latex_output()) badger::badge_doi("10.5281/zenodo.5083309", "blue")

Installation

You can(not) install the released version of codos from CRAN with:

install.packages("codos")

And the development version from GitHub with:

# install.packages("devtools")
devtools::install_github("special-uor/codos", "dev")

cat(paste0("
-----

Note: Some of the equations on this document are not displayed properly (due to a server issue), check out the [README-extended.pdf](README-extended.pdf).

-----
"))

Background:

Vapour-pressure deficit (`vpd`)

vpd is given by mean daily growing season temperature, tmp [°C] and moisture index, mi [-]. Using the CRU TS 4.04 dataset [@cru404] we found the following relation:

$$\text{vpd} = 4.6123 \times \exp(0.0609 \times \text{tmp}-0.8726 \times \text{mi})$$

The steps performed were:

Generate a monthly climatology for the period between 1961 and 1990 (inclusive). Variables used: cld, pre, tmn, tmx, vap.

r # Monthly climatology for `tmn` codos::monthly_clim("cru_ts4.04.1901.2019.tmn.dat.nc", "tmn", 1961, 1990)

Output file:

bash "cru_ts4.04.1901.2019.tmn.dat-clim-1961-1990.nc"

Interpolate the monthly data to daily. Variables used: cld, pre, tmn, tmx, vap.

r # Monthly to daily interpolation for `tmn` codos::nc_int("cru_ts4.04.1901.2019.tmn.dat-clim-1961-1990.nc", "tmn")

Output file:

bash "cru_ts4.04.1901.2019.tmn.dat-clim-1961-1990-int.nc"

Calculate daily temperature, tmp. Variables used: tmn and tmx.

r codos::daily_temp(tmin = list(filename = "cru_ts4.04.1901.2019.tmn.dat-clim-1961-1990-int.nc", id = "tmn"), tmax = list(filename = "cru_ts4.04.1901.2019.tmx.dat-clim-1961-1990-int.nc", id = "tmx"), output_filename = "cru_ts4.04-clim-1961-1990-daily.tmp.nc")

Calculate mean growing season for daily temperature

r codos::nc_gs("cru_ts4.04-clim-1961-1990-daily.tmp.nc", "tmp", thr = 0)

Output file:

bash "cru_ts4.04-clim-1961-1990-daily.tmp-gs.nc"

Calculate potential evapotranspiration (pet)

Install SPLASH (unofficial R package) as follows:

r remotes::install_github("villegar/splash", "dev")

Or, download from the official source: https://bitbucket.org/labprentice/splash.

```r elv <- codos:::nc_var_get("halfdeg.elv.nc", "elv")$data tmp <- codos:::nc_var_get("cru_ts4.04.1901.2019.daily.tmp.nc", "tmp")$data cld <- codos:::nc_var_get("cru_ts4.04.1901.2019.cld.dat-clim-1961-1990-int.nc", "cld")$data

codos::splash_evap(output_filename = "cru_ts4.04-clim-1961-1990-pet.nc", elv, # Elevation, 720x360 grid sf = 1 - cld / 100, tmp, year = 1961, # Reference year lat = codos::lat, lon = codos::lon) ```

Output file:

bash "cru_ts4.04-clim-1961-1990-pet.nc"

Calculate moisture index (mi)

$$MI_{i,j} = \frac{\text{Total precipitation}}{\text{Total PET}}$$

r pet <- codos:::nc_var_get("cru_ts4.04-clim-1961-1990-pet.nc", "pet")$data pre <- codos:::nc_var_get("cru_ts4.04.1901.2019.pre.dat-new-clim-1961-1990-int.nc", "pre")$data codos::nc_mi(output_filename = "cru_ts4.04-clim-1961-1990-mi.nc", pet, # potential evapotranspiration pre) # precipitation

Output file:

bash "cru_ts4.04-clim-1961-1990-mi.nc"

Approximate vpd

r tmp <- codos:::nc_var_get("cru_ts4.04-clim-1961-1990-daily.tmp.nc", "tmp")$data vap <- codos:::nc_var_get("cru_ts4.04.1901.2019.vap.dat-clim-1961-1990-int.nc", "vap")$data output_filename <- file.path(path, "cru_ts4.04-clim-1961-1990-vpd-tmp.nc") codos::nc_vpd(output_filename, tmp, vap)

Output file:

bash "cru_ts4.04-clim-1961-1990-vpd-tmp.nc"

Find the coeffients for the following equation

$$\text{vpd} = a \times \exp(\text{kTmp} \times \text{tmp}-\text{kMI} \times \text{mi})$$

mi <- codos:::nc_var_get("cru_ts4.04-clim-1961-1990-mi.nc", "mi")$data
Tmp <- codos:::nc_var_get("cru_ts4.04-clim-1961-1990-daily.tmp-gs.nc", "tmp")$data
vpd <- codos:::nc_var_get("cru_ts4.04-clim-1961-1990-vpd-tmp-gs.nc", "vpd")$data

# Apply ice mask
mi[codos:::ice_mask] <- NA
Tmp[codos:::ice_mask] <- NA
vpd[codos:::ice_mask] <- NA

# Filter low temperatures, Tmp < 5
mi[Tmp < 5] <- NA
Tmp[Tmp < 5] <- NA

# Create data frame
df <- tibble::tibble(Tmp = c(Tmp), 
                     vpd = c(vpd),
                     MI = c(mi))

# Filter grid cells with missing Tmp, vpd, or MI
df <- df[!is.na(df$Tmp) & !is.na(df$vpd) & !is.na(df$MI), ]

# Linear approximation
lmod <- lm(log(vpd) ~ Tmp + MI, data = df)
# Non-linear model
exp_mod <- nls(vpd ~ a * exp(kTmp * Tmp - kMI * MI),
               df,
               start = list(a = exp(coef(lmod)[1]),
                            kTmp = coef(lmod)[2],
                            kMI = coef(lmod)[3]),
               control = list(maxiter = 200))

Summary statistics:

df <- read.csv("documentation/codos/inst/extdata/df.csv")
# Linear approximation
lmod <- lm(log(vpd) ~ Tmp + MI, data = df)
# Non-linear model
exp_mod <- nls(vpd ~ a * exp(kTmp * Tmp - kMI * MI),
               df,
               start = list(a = exp(coef(lmod)[1]),
                            kTmp = coef(lmod)[2],
                            kMI = coef(lmod)[3]),
               control = list(maxiter = 200))

summary(exp_mod)
coefficients(exp_mod)

knitr::asis_output('\n\n\\newpage')

Corrected `mi` from reconstructed `mi`

The following equations were used:

$$f = e/1.6 = D/[\text{c}_\text{a}(1-\chi)]$$

$$\chi = \frac{\xi \times \text{vpd}^{1/2} + \text{vpd}}{\text{c}\text{a}/\left(\text{c}\text{a}+9.7\right)}$$

$$\xi = [\beta (K + \Gamma^) / (1.6 \eta^)] ^ {1/2}$$

where:

$e$ ratio of water lost to carbon fixed [–]
$\text{vpd}$ vapour pressure deficit [Pa]
$\text{c}_\text{a}$ ambient CO2 partial pressure [Pa]
$\chi$ ratio of leaf-internal to ambient CO2 partial pressures [–]
$\xi$ stomatal sensitivity factor [Pa1/2]
$\Gamma^*$ photorespiratory compensation point [Pa]: a function of temperature and elevation
$\beta$ ratio of cost factors for carboxylation and transpiration = 146 [–]
$K$ effective Michaelis constant of Rubisco [Pa]: a function of temperature and elevation
$\eta^*$ viscosity of water relative to its value at 25°C [–]

And the equilibrium relation:

$$f(\text{T}\text{c1}, \text{MI}\text{1}, \text{c}\text{a,1}) = f(\text{T}\text{c0}, \text{MI}\text{0}, \text{c}\text{a,0})$$ where:

$\text{T}_\text{c1}$ past temperature (assume equal to reconstructed value) [K]
$\text{MI}_\text{1}$ past MI (unknown) [–]
$\text{c}_\text{a,1}$ past ambient CO2 partial pressure [Pa], adjusted for elevation
$\text{T}_\text{c0}$ present temperature [K]
$\text{MI}_\text{0}$ reconstructed MI [–]
$\text{c}_\text{a,1}$ 'recent' ambient CO2 partial pressure [Pa], adjusted for elevation

Steps in the solution:

Evaluate $f(\text{T}\text{c0}, \text{MI}\text{0}, \text{c}_\text{a,0})$
Equate this to:

$$[\xi(\text{T}\text{c1}, z) \times \text{vpd}_1 ^ {1/2} + \text{vpd}_1] / [\text{c}\text{a,1}(z) / (\text{c}_\text{a,1}(z) + 9.7)] $$

where:

$z$ is elevation
$\text{vpd}_1$ is past vapour pressure deficit

And solve for $\text{vpd}_1$.

Convert $\text{vpd}1$ back to MI (at temperature $\text{T}\text{c1}$), to yield an estimate of $\text{MI}_\text{1}$.

Using codos, all the steps translate to a simple function call

corrected_mi <- codos::corrected_mi(present_t,
                                    past_temp,
                                    recon_mi,
                                    modern_co2,
                                    past_co2)

Note that this function takes temperatures in [°C] and ambient CO$_2$ partial pressures in [$\mu\text{mol}/\text{mol}$] (unless, scale_factor is overwritten, e.g. scale_factor = 1 to use ambient CO$_2$ partial pressures in [Pa]).

More details: