geco-bern / rpmodel

An R implementation of the P-model.

Purpose

rpmodel provides an implementation of the P-model (Prentice et al., 2014; Wang et al., 2017; Stocker et al., 2020) for predicting acclimated photosynthetic parameters, assimilation, and dark respiration rates as a function of the environment. The main function is rpmodel() which returns a list of variables that are mutually consistent within the theory of the P-model (see Usage ). Further functions used within rpmodel() are also provided through the package.

Important note:

The P-model predicts how photosynthesis acclimates to a changing environment, coordinating stomatal conductance, Vcmax and Jmax. This yields a model that has the form of a light use efficiency model, where gross primary production scales linearly with absorbed light, as described in Stocker et al. 2020. It is important to note that this implies that the P-model is valid only for simulating responses to the environment that evolve over the time scale at which the photosynthetic machinery (e.g., Rubisco) can be assumed to acclimate. Sensible choices are on the order of a couple of weeks to a month. In other words, the arguments (climatic forcing), provided to rpmodel() should represent typical daytime mean values, averaged across a couple of weeks. The output is then representative also for average values across the same time scale.

Usage

This loads the rpmodel package and executes the rpmodel() function without $J_{\text{max}}$ limitation (argument method_jmaxlim = "none"), and with a temperature-independent quantum yield efficiency (argument do_ftemp_kphio = FALSE):

library(rpmodel)
out_pmodel <- rpmodel( 
  tc             = 20,           # temperature, deg C
  vpd            = 1000,         # Pa,
  co2            = 400,          # ppm,
  fapar          = 1,            # fraction  ,
  ppfd           = 300,          # mol/m2/d,
  elv            = 0,            # m.a.s.l.,
  kphio          = 0.049977,     # quantum yield efficiency as calibrated for setup ORG by Stocker et al. 2020 GMD,
  beta           = 146,          # unit cost ratio a/b,
  c4             = FALSE,
  method_optci   = "prentice14",
  method_jmaxlim = "none",
  do_ftemp_kphio = FALSE,        # corresponding to setup ORG
  do_soilmstress = FALSE,        # corresponding to setup ORG
  verbose        = TRUE
  )

For more information and examples see Usage.

Installation

Stable release

rpmodel is available on CRAN here. To install and load, run the following commands in your R terminal:

install.packages("rpmodel")
library(rpmodel)

Development release

To install and load the latest version of the rpmodel package (development release, not yet on CRAN) run the following command in your R terminal:

if(!require(devtools)){install.packages(devtools)}
devtools::install_github( "stineb/rpmodel", build_vignettes = TRUE )
library(rpmodel)

Author and contact

Benjamin Stocker [email protected]

References

Stocker, B. D., Wang, H., Smith, N. G., Harrison, S. P., Keenan, T. F., Sandoval, D., Davis, T., and Prentice, I. C.: P-model v1.0: an optimality-based light use efficiency model for simulating ecosystem gross primary production, Geosci. Model Dev., 13, 1545–1581, https://doi.org/10.5194/gmd-13-1545-2020, 2020.

Wang, H., Prentice, I. C., Keenan, T. F., Davis, T. W., Wright, I. J., Cornwell, W. K.,Evans, B. J., and Peng, C.: Towards a universal model for carbon dioxide uptake by plants, Nat Plants, 3, 734–741, 2017.

Prentice, I. C., Dong, N., Gleason, S. M., Maire, V., and Wright, I. J.: Balancingthe costs of carbon gain and water transport: testing a new theoretical frameworkfor plant functional ecology, Ecology Letters, 17, 82–91, 10.1111/ele.12211, 2014.

Acknowledgement

This project was funded by Marie Sklodowska-Curie fellowship H2020-MSCA-IF-2015, project FIBER, grant number 701329.

rpmodel's Issues

Numerical instability in density_h20

Sorry about these bugreps! I'm chasing down some numerical warnings in my python implementation that cropped up using real world inputs with more extreme ranges in test inputs.

There is particular behaviour in the equations for density_h20 that create aberrant results. Specifically, temperatures between about -44 and -46 show extreme numerical instability. The code below shows how the mean and sd of the predicted density go haywire for small changes in temperature in this region:

t <- seq(-88, 58, length=1001)
r <- runif(100, -0.5,0.5)
sd_t <- sapply(t, function(xx) sd(density_h2o(xx + r, p=101325)))
plot(sd_t ~ t, type='l')
rho_t <- sapply(t, function(xx) mean(density_h2o(xx + r, p=101325)))
plot(rho_t ~ t, type='l')

This isn't really a flaw in the equation - it is just being extrapolated way outside the intended range of use and happens to have behaviour that causes all sorts of downstream numerical issues.

I think the best fix here is probably a more general constraint on temperature inputs to the pmodel - so that no predictions are made for sub-zero temperatures, given how little the freezing point changes under environmental temperature and pressures.

Re-submit to CRAN

prepare for release:

• Check CRAN issues (ISSUES)
• Update DESCRIPTION - version increment!
• Update / check NEWS.md
• Update cran-comments.md for release
• urlchecker::url_check()
• devtools::check(document = FALSE)
• rhub::rhub_setup() [if not done yet, may have to do rhub::rhub_doctor(), too]
• Check build:

# Build package
pkg <- devtools::build()

# Check package built from source
devtools::check_built(path = pkg)

• R CMD check --as-cran

system("R CMD check --as-cran ~/rpmodel_1.2.3.tar.gz")

• git push remaining changes

Submit to CRAN:

• devtools::submit_cran()
• Approve email
• Accepted email
• Update release

patm calculation may better use kT0 = 288.15K (see Berberan-Santos et al., 1997; Equ.40)

Dear Beni,

According to the original elevation dependence of atmospheric pressure formula proposed by (Berberan-Santos et al., 1997; Equ.40), kT0 better to be set as 288.15K rather than 298.15K, since it better represent the global annual mean surface temperature of the Earth. 288.15K was also originally proposed by Berberan-Santos et al., 1997, in its Equation 40.

To solve this, it is suggested to change kTo <- 288.15K, in patm function at rpmodel/R/subroutines.R.

Where parallel to this issue, also changing kTo <- 288.15K in calc_patm function at ingestr/R/calc_patm.R, since we will rely on this function to download vpd in site-simulations.

Much thanks
Yunke

Inequality in a_j == a_c calculation when iabs == 0

Hi,

Connected to #11, one way you get a mismatch between a_j and a_c is when iabs is zero.

a_j

fact_jmaxlim <- vcmax * (ci + 2 * gammastar)/(kphio * iabs * (ci + kmm))
a_j <- kphio * iabs * (ci - gammastar)/(ci + 2 * gammastar) * fact_jmaxlim

Because fact_jmaxlim has vcmax in the numerator and iabs in the denominator, this gives 0/0 and hence both it and a_j are NaN when iabs is zero.

a_c

In contrast:

a_c <- vcmax * (ci - gammastar)/(ci + kmm)

Here iabs only features in numerator via vcmax and so is 0/real = 0.

In practice, I guess the equation for a_j should be set to zero when there is no photosynthetically active light, and that needs to be trapped to allow the test in #11 to work correctly.

Atmospheric CO2 driver not vectorized

Unlike other inputs, if a sequence of CO2 values is provided, only the first value is used.

Expected behavior (this is what currently happens with, e.g., a VPD sweep):

> rpmodel( tc = 25, vpd = seq(500, 3000, 500), co2 = 400, fapar = 1, ppfd = 10, elv = 0, do_soilmstress = FALSE, do_ftemp_kphio = FALSE)$gpp
[1] 2.224166 2.115682 2.035845 1.970669 1.914766 1.865385

Actual behavior (when a CO2 sweep is used instead):

> rpmodel( tc = 25, vpd = 500, co2 = seq(350, 400, 10), fapar = 1, ppfd = 10, elv = 0, do_soilmstress = FALSE, do_ftemp_kphio = FALSE)$gpp
[1] 2.027605

Wrong output size

Using vectors for all variables (temperature, vpd, co2, ppfd, fapar, patm) retuned the wrong number of outputs (something like 2 millions lines even though input vectors had a size of 1460).
A modification was made to line 336 of rpmodel.R . It now reads:
iabs <- fapar * ppfd
I don't know if it would work for all cases.

Stomatal conductance in C4 plants

Stomatal conductance is currently calculated using the dummy value 1.0 for optimal chi. For calculating g_s, the true internal CO2 partial pressure is needed. This can be calculated for C4 plants using the same equations but a different beta cost, which Colin has suggested should be:

Both Lin et al 2015 and DeKauwe et al 2015 provide estimates for the $g_1$
parameter for C3 and C4 plants, with a ratio of C3/C4 values of around 3. The
$g_1$ parameter is equivalent to $\xi$ in the P model. Given that
$\xi \propto \surd\beta$, a reasonable default for C4 plants is that
$\beta = 146 / 9 \approx 16.222$

So, a new optimal chi function is needed that returns $\chi$ for use in appropriate calculations.

Does the latest Pmodel support grid data?

Does the latest Pmodel support grid data?

Bug in do_ftemp_kphio

The same issue as in #16 also affect the mechanics of do_ftemp_kphio:

https://github.com/computationales/rpmodel/blob/master/R/rpmodel.R#L287

  kphio <- ifelse(
    do_ftemp_kphio,
    ftemp_kphio( tc, c4 ) * kphio,
    kphio
  )

The argument do_ftemp_kphio is a logical value, so this structure guarantees that kphio will be truncated when do_ftemp_kphio is TRUE:

> kphio <- ifelse(TRUE, c('lots', 'of', 'numbers'), c('some', 'other','numbers'))                                                                             
> kphio                                                                                                                                               
[1] "lots"

I am working on fixes for these bugs - I'll submit some pull requests!

Bug in C4 calculation

There is a problem with the test for assimilation == GPP when using C4=TRUE:

> ret <- rpmodel(tc = 20, vpd = 1000, co2 = 400, ppfd = 30, elv = 0, c4=TRUE, fapar=1)
Error in rpmodel(tc = 20, vpd = 1000, co2 = 400, ppfd = 30, elv = 0, c4 = TRUE,  : 
  rpmodel(): Assimilation and GPP are not identical.
> ret <- rpmodel(tc = 20, vpd = 1000, co2 = 400, ppfd = 30, elv = 0, c4=FALSE, fapar=1)
>

The lines are:

    assim <- ifelse(a_j < a_c, a_j, a_c)
    if (any(abs(assim - gpp/c_molmass) > 0.001)) 
        stop("rpmodel(): Assimilation and GPP are not identical.")

I can't immediately see what is wrong - it isn't NA issues as in #11 and the values genuinely are unequal. Having said that, I think all.equal is more robust here as well.

Also - the first line there creating assim: the previous code has just tested that a_j == a_c, so presumably a_j or a_c could be used directly rather than going to the computational expense of combining them?

Bug in soilmstress

Hi,

There is a bug in soilmstress at the line:
https://github.com/computationales/rpmodel/blob/master/R/subroutines.R#L124

The problem only arises if a user provides soilm as a scalar but meanalpha as an array. That's an odd thing to want to do, but might arise if someone wants to see how the stress factor changes for a given soil moisture in different sites.

In this case - in the nested ifelse statements - the first test 'iterates' over the comparison of the single scalar soilm to x1. If meanalpha is an array, the output is therefore either 1 or the first element only of the no argument:

Browse[2]> ifelse(test = 2 > 4, yes = 1, no = runif(10))                                                                                                                         
[1] 0.5927906

A bit obscure, I admit!

Forced scaling in do_ftemp_kphio with C4 plants

Hi,

There is something odd going on with the code for setting kphio here:
https://github.com/computationales/rpmodel/blob/master/R/rpmodel.R#L305

  if (do_ftemp_kphio) {
    kphio <- ftemp_kphio( tc, c4 ) * kphio
  } else {
    if (c4){
      kphio <- ftemp_kphio( 15.0, c4 ) * kphio
    }
  }

If the temperature correction is turned off, C4 plants are forced to be reduced by a scaling factor ftemp_kphio( 15.0, c4 ) = 0.2816. That seems like an error?

NA handling bug in testing a_j = a_c

Hi,

This line introduces a problem with NA values:

https://github.com/stineb/rpmodel/blob/master/R/rpmodel.R#L437

In the presence of NAs (but no other problems) the test fails:

> any(c(FALSE, NA))  # returns NA
[1] NA
> if(any(c(FALSE, NA))) { stop('a_j!=a_c')}  # if() statement chokes,
Error in if (any(c(FALSE, NA))) { : missing value where TRUE/FALSE needed

With grid maps as inputs, NAs are pretty likely and are handled by all steps up to this one. One alternative here is all.equal, which handles NAs and allows a tolerance to be set:

> all.equal(c(1,2,3,4, NA), c(1,2,3,4.001, NA), tol=0.001)
[1] TRUE
> all.equal(c(1,2,3,4, NA), c(NA, 2,3,4.001, NA), tol=0.001)
[1] "'is.NA' value mismatch: 2 in current 1 in target"

It would probably be useful to handle $a_j != a_c$ and mismatched NA values separately though? So trap the output of that all.equal and report separately or perhaps use a separate xor(is.na(a_j), is.na(a_c))

Negative value corrections in calc_ftemp_kphio

Hi Beni,

There's an issue with the correction to filter out negative values in calc_ftemp_kphio in version 1.1.2 . The max() function finds the global maximum of the inputs, which of course collapses matrix and array inputs down to a scalar (and gets hung up on missing values and NaN).

> tc <- array(c(NA, -9:28, NaN), dim=c(4,5,2))
> calc_ftemp_kphio(tc)
[1] NA

Substituting pmax(ftemp, 0) preserves the shape of ftemp and calculates maxima for each matrix/array entry.

Scaling of do_ftemp_kphio for C4 plants

The values used here:

https://github.com/computationales/rpmodel/blob/master/R/subroutines.R#L358

ftemp = -0.008 + 0.00375 * tc - 0.58e-4 * tc^2

are incorrect. Shirley has pointed out that the values have been scaled twice in error and that these should be:

ftemp = -0.064 + 0.03 * tc - 0.000464 * tc^2 

correct calc_optimal_chi for dew conditions (vpd<0)

Hello Beni,
I found a condition where your model crashes. On fluxnet sites, I compute VPD from Qair, press and Temperature.
In dew conditions we can have negative vpd, with qair>qsat.
In calc_optimal_chi you have sqrt(vpd), which leads to NaN values for vpd<0. You can easily correct that by using sqrt(max(vpd,0)).
Cheers,

Stocker2020 Eq.12

Dear Stocker,

Following your 2020 "P-model v1.0: an optimality-based light use efficiency model for simulating ecosystem gross primary production", the deduced Eq. 12 different from yours.
Could you have a look please?

The following is the corresponding MATLAB code:

syms gamma c_i c_s C_i C_a A g_sw beta Gamma R_d Chi g_0 epsilon xi eta K D
% epsilo = 1.6

% K = 2 * Gamma
xi = sqrt( beta * (K + Gamma) / (epsilon * eta ))
Chi = Gamma/C_a + (1 - Gamma/C_a) / xi/(xi + sqrt(D))

C_i = simplify( Chi * C_a )

m = (C_i - Gamma) / (C_i + 2 * Gamma);
m = simplify(m)