BioCro II Paper: Section 1.1 Example

Introduction

This is a demonstration of the example discussed in Section 1.1 and Appendix 1 of the BioCro II paper (Lochocki et al. 2022). A few corrections to and clarifications of the exposition given in that paper are given in the comments below.

The Code

library(BioCro)
library(lattice)
library(knitr) # for kable(), which yields nicer-looking tables

## customized version of kable with defined options:
columns_to_print <- c('time', 'Q', 'mass_gain', 'Root', 'Leaf')
show_row_number <- FALSE
format <- list(scientific = FALSE, digits = 7)
cable <- function(x, ...) {
   kable(x[, columns_to_print],
         row.names = show_row_number,
         format.args = format,
         ...)
}

## Set plotting character globally
trellis.par.set("plot.symbol", list(pch = '.'))

################################################################################

parameters <- list(
    alpha_rue = 0.07, # kg / mol
    SLA = 25,         # m^2 / kg
    C_conv = 0.03,    # kg / mol
    f_leaf = 0.2,     # kg / kg
    f_root = 0.8,     # kg / kg
    timestep = 1      # s
)

initial_values <- list(
    Leaf = 1, # kg
    Root = 1  # kg
)

Q <- function(time) sin(time/3600/12 * pi) * 2000e-6 # mol / m^2 / s

times <- 0:(3600 * 12) # seconds

light_intensity <- data.frame(
    time = times,
    Q = Q(times)
)

result <- run_biocro(
    initial_values,
    parameters,
    light_intensity,
    'BioCro:example_model_mass_gain',
    'BioCro:example_model_partitioning'
)

xyplot(Q ~ time, data = result, xlab = "time (s)", ylab = "Q (mol / m^2 / s)")

xyplot(Leaf ~ time, data = result, xlab = "time (s)", ylab = "Leaf (kg)")

A modification: 24 hour simulation

Although in the paper it is stated that “Light intensity is given as a table of values at every second in a 24-h period,”, the period here is actually only 12 hours (time ranges from 0 to 3600 × 12 seconds = 12 hours). As noted above, for the purposes of this example, it may be assumed that time ranges from sunrise to sunset of a 12-hour day so that “time” represents the number of seconds that have elapsed since sunrise.

Alternatively, we could have assumed time to range over a 24-hour period starting and ending at solar midnight, with Q being zero before sunrise and after sunset. This entails a slight modification to the defining equations:

times <- 0:(3600 * 24) # seconds

Q <- function(time) pmax(0, sin((time/3600 - 6)/12 * pi) * 2000e-6) # mol / m^2 / s

(The “− 6” here adjusts the phase of the sine function so that the positive values correspond to the daylight hours, assumed to begin at 6 a.m. solar time.)

Here is how the 24-hour version of the example would look. A few rows of the result table are printed at significant portions of the day: at the beginning, just after midnight, when it is dark and there is no mass accumulation; at dawn, when the rate of mass accumulation first becomes non-zero; at noon, when the flux Q is highest; at dusk, when the rate of mass accumulation again drops to zero; and at the end of the simulation, when the final leaf and root mass are known.

Note that the tabular result showed below omits some columns from the result, such as ncalls, which is not of particular interest here. And the values in the time column should be interpreted as the number of elapsed seconds from solar midnight (not the number of (fractional) hours that have elapsed since midnight at the start of January 1, which is the correct interpretation in a “normal” BioCro simulation).

times <- 0:(3600 * 24) # seconds

Q <- function(time) {
    pmax(0, sin((time/3600 - 6)/12 * pi) * 2000e-6) # mol / m^2 / s
}

light_intensity <- data.frame(
    time = times,
    Q = Q(times)
)

result <- run_biocro(
    initial_values,
    parameters,
    light_intensity,
    'BioCro:example_model_mass_gain',
    'BioCro:example_model_partitioning'
)

## cable = customized version of kable (see above)
## kable = (nicely-formatted) knitr table
cable(result[1:4,]) # beginning

cable(result[seq(21591, 21661, 10),]) # around dawn

cable(result[43199:43203,]) # around mid day

cable(result[seq(64741, 64811, 10),]) # around dusk

cable(result[86398:86401,]) # end

In the following graphs, we have changed the labeling on the x axis so that time is shown in hours rather than seconds.

xyplot(Q ~ time / 3600, data = result, xlim = c(0, 24), xlab = "time (in hours)", ylab = "Q (mol / m^2 / s)")

xyplot(Leaf ~ time / 3600, data = result, xlim = c(0, 24), xlab = "time (in hours)", ylab = "Leaf (kg)")

Modules with time as an input

The usual assumption in BioCro’s crop growth models is that the time quantity represents the time, in hours, that has elapsed from the beginning of the year, not the number of seconds that have elapsed from the start of the simulation.

It follows, then, that the way we have been using the time variable in the previous examples in inconsistent with the usual interpretation of time in two respects: First, the units were in seconds rather than days. Second, the time represented the amount of time elapsed from the beginning of the simulation rather than an “absolute” time, that is, a time value specifying a specific time of day on a specific day of the year (a date and time on the calendar).

Thus, direct modules that rely on the usual meaning for time, such as BioCro:format_time, could not be used as part of the system considered in the examples above.

Yet, there is no issue with using the time variable in this unconventional way, as long as we know what we are doing. In fact, other BioCro module libraries and models are free to use any time convention.