Package design
PlantBiophysics.jl
is designed to ease the computations of biophysical processes in plants and other objects. It uses PlantSimEngine.jl
, so it shares the same ontology (same concepts and terms).
Definitions
Processes
A process is defined in PlantSimEngine as a biological or a physical phenomena. At this time PlantBiophysics.jl
implements four processes:
- light interception
- energy balance
- photosynthesis
- stomatal conductance
Models
A process is simulated using a particular implementation of a model. Each model is implemented using a structure that lists the parameters of the model. For example, PlantBiophysics provides the Beer
structure for the implementation of the Beer-Lambert law of light extinction.
You can see the list of available models for each process in the sections about the models, e.g. here for photosynthesis.
Models can use three types of entries:
- Parameters
- Meteorological information
- Variables
Parameters are constant values that are used by the model to compute its outputs. Meteorological information are values that are provided by the user and are used as inputs to the model. Variables are computed by the model and can optionally be initialized before the simulation.
Users can choose which model is used to simulate a process using the ModelList
structure from PlantSimEngine. ModelList
is also used to store the values of the parameters, and to initialize variables.
Let's instantiate a ModelList
with the Beer-Lambert model of light extinction. The model is implemented with the Beer
structure and has only one parameter: the extinction coefficient (k
).
using PlantSimEngine, PlantBiophysics
ModelList(Beer(0.5))
PlantSimEngine.DependencyGraph{Dict{Symbol, PlantSimEngine.SoftDependencyNode}}(Dict{Symbol, PlantSimEngine.SoftDependencyNode}(:light_interception => Beer{Float64}
), Dict{Symbol, DataType}())TimeStepTable{Status{(:LAI, :aPPFD), Tuple{...}(1 x 2):
╭─────┬─────────┬─────────╮
│ Row │ LAI │ aPPFD │
│ │ Float64 │ Float64 │
├─────┼─────────┼─────────┤
│ 1 │ -Inf │ -Inf │
╰─────┴─────────┴─────────╯
What happened here? We provided an instance of a model to a ModelList
that automatically associates it to the process it simulates (i.e. the light interception).
We see that we only instantiated the ModelList
for the light extinction process. What about the others like photosynthesis or energy balance ? Well there is no need to give models if we have no intention to simulate them.
Parameters
A parameter is a constant value that is used by a model to compute its outputs. For example, the Beer-Lambert model uses the extinction coefficient (k
) to compute the light extinction. The Beer-Lambert model is implemented with the Beer
structure, which has only one field: k
. We can see that using fieldnames
:
fieldnames(Beer)
(:k,)
Some models are shipped with default values for their parameters. For example, the Monteith
model that simulates the energy balance has a default value for all its parameters. Here are the parameter names:
fieldnames(Monteith)
(:aₛₕ, :aₛᵥ, :ε, :maxiter, :ΔT)
And their default values:
Monteith()
Monteith{Float64, Int64}(2, 1, 0.955, 10, 0.01)
But if we need to change the values of some parameters, we can usually give them as keyword arguments:
Monteith(maxiter = 100, ΔT = 0.001)
Monteith{Float64, Int64}(2, 1, 0.955, 100, 0.001)
Perfect! Now is that all we need to make a simulation? Well, usually no. Models need parameters, but also input variables.
Variables (inputs, outputs)
Variables are computed by models, and can optionally be initialized before the simulation. Variables and their values are stored in the ModelList
, and are initialized automatically or manually.
Hence, ModelList
objects store two fields:
fieldnames(ModelList)
(:models, :status, :vars_not_propagated)
The first field is a list of models associated to the processes they simulate. The second, :status
, is used to hold all inputs and outputs of our models, called variables. For example the Beer
model needs the leaf area index (LAI
, m^{2} \cdot m^{-2}) to run.
We can see which variables are needed as inputs using inputs
from PlantSimEngine:
using PlantSimEngine
inputs(Beer(0.5))
(:LAI,)
We can also see the outputs of the model using outputs
from PlantSimEngine:
using PlantSimEngine
outputs(Beer(0.5))
(:aPPFD,)
If we instantiate a ModelList
with the Beer-Lambert model, we can see that the :status
field has two variables: LAI
and PPDF
. The first is an input, the second an output.
using PlantSimEngine, PlantBiophysics
m = ModelList(Beer(0.5))
keys(m.status)
(:LAI, :aPPFD)
To know which variables should be initialized, we can use to_initialize
from PlantSimEngine:
m = ModelList(Beer(0.5))
to_initialize(m)
(light_interception = (:LAI,),)
Their values are uninitialized though (hence the warnings):
(m[:LAI], m[:aPPFD])
([-Inf], [-Inf])
Uninitialized variables have often the value returned by typemin()
, e.g. -Inf
for Float64
:
typemin(Float64)
-Inf
Prefer using to_initialize
rather than inputs
to check which variables should be initialized. inputs
returns the variables that are needed by the model to run, but to_initialize
returns the variables that are needed by the model to run and that are not initialized. Also to_initialize
is more clever when coupling models (see below).
We can initialize the variables by providing their values to the status at instantiation:
m = ModelList(Beer(0.5), status = (LAI = 2.0,))
PlantSimEngine.DependencyGraph{Dict{Symbol, PlantSimEngine.SoftDependencyNode}}(Dict{Symbol, PlantSimEngine.SoftDependencyNode}(:light_interception => Beer{Float64}
), Dict{Symbol, DataType}())TimeStepTable{Status{(:LAI, :aPPFD), Tuple{...}(1 x 2):
╭─────┬─────────┬─────────╮
│ Row │ LAI │ aPPFD │
│ │ Float64 │ Float64 │
├─────┼─────────┼─────────┤
│ 1 │ 2.0 │ -Inf │
╰─────┴─────────┴─────────╯
Or after instantiation using init_status!
(from PlantSimEngine):
m = ModelList(Beer(0.5))
init_status!(m, LAI = 2.0)
[ Info: Some variables must be initialized before simulation: (light_interception = (:LAI,),) (see `to_initialize()`)
We can check if a component is correctly initialized using is_initialized
(from PlantSimEngine):
is_initialized(m)
true
Some variables are inputs of models, but outputs of other models. When we couple models, we have to be careful to initialize only the variables that are not computed.
Climate forcing
To make a simulation, we usually need the climatic/meteorological conditions measured close to the object or component.
The PlantMeteo.jl
package provides a data structure to declare those conditions, and to pre-compute other required variables. The most basic data structure is a type called Atmosphere
, which defines the conditions for a steady-state, i.e. the conditions are considered at equilibrium. Another structure is available to define different consecutive time-steps: Weather
.
The mandatory variables to provide for an Atmosphere
are: T
(air temperature in °C), Rh
(relative humidity, 0-1), Wind
(the wind speed in m s-1) and P
(the air pressure in kPa). We can declare such conditions like so:
using PlantMeteo
meteo = Atmosphere(T = 20.0, Wind = 1.0, P = 101.3, Rh = 0.65)
Atmosphere(date = Dates.DateTime("2024-07-11T16:04:58.017"), duration = Dates.Second(1), T = 20.0, Wind = 1.0, P = 101.3, Rh = 0.65, Precipitations = 0.0, Cₐ = 400.0, e = 1.5255470730405223, eₛ = 2.3469954969854188, VPD = 0.8214484239448965, ρ = 1.2037851579511918, λ = 2.4537e6, γ = 0.06723680111943287, ε = 0.5848056484857892, Δ = 0.14573378083416522, clearness = Inf, Ri_SW_f = Inf, Ri_PAR_f = Inf, Ri_NIR_f = Inf, Ri_TIR_f = Inf, Ri_custom_f = Inf)
More details are available from the dedicated section.
Simulation
Simulation of processes
Making a simulation is rather simple, we simply use the run!
function provided by PlantSimEngine
:
run!(model_list, meteo)
The first argument is the model list (see ModelList
from PlantSimEngine
), and the second defines the micro-climatic conditions (more details below in Climate forcing).
The ModelList
should be initialized for the given process before calling run!
. See Variables (inputs, outputs) for more details.
Example simulation
For example we can simulate the stomatal_conductance
of a leaf like so:
using PlantMeteo, PlantSimEngine, PlantBiophysics
meteo = Atmosphere(T = 20.0, Wind = 1.0, P = 101.3, Rh = 0.65)
leaf = ModelList(
Medlyn(0.03, 12.0),
status = (A = 20.0, Dₗ = meteo.VPD, Cₛ = 400.0)
)
run!(leaf, meteo)
leaf[:Gₛ]
1-element Vector{Float64}:
0.7420047415309556
Outputs
The status
field of a ModelList
is used to initialize the variables before simulation and then to keep track of their values during and after the simulation. We can extract the simulation outputs of a model list using the status
function (from PlantSimEngine).
The status can either be a Status
type if simulating only one time-step, or a TimeStepTable
(from PlantMeteo
) if several.
Let's look at the status of our previous simulated leaf:
We can extract the value of one variable using the status
function, e.g. for the stomatal conductance:
status(leaf, :Gₛ)
1-element Vector{Float64}:
0.7420047415309556
Or similarly using the dot syntax:
leaf.status.Gₛ
1-element Vector{Float64}:
0.7420047415309556
Or much simpler (and recommended), by indexing directly the model list:
leaf[:Gₛ]
1-element Vector{Float64}:
0.7420047415309556
Another simple way to get the results is to transform the outputs into a DataFrame
:
using DataFrames
DataFrame(leaf)
Row | Dₗ | Cₛ | A | Gₛ | timestep |
---|---|---|---|---|---|
Float64 | Float64 | Float64 | Float64 | Int64 | |
1 | 0.821448 | 400.0 | 20.0 | 0.742005 | 1 |
The output from DataFrame
is adapted to the kind of simulation you did: one row per time-steps, and per component models if you simulated several.
Model coupling
A model can work either independently or in conjunction with other models. For example a stomatal conductance model is often associated with a photosynthesis model, i.e. it is called from the photosynthesis model.
Several models proposed in PlantBiophysics.jl
are hard-coupled models, i.e. one model calls another. For example, the Fvcb
structure is the implementation of the Farquhar–von Caemmerer–Berry model for C3 photosynthesis (Farquhar et al., 1980; von Caemmerer and Farquhar, 1981) calls a stomatal conductance model. Hence, using Fvcb
requires a stomatal conductance model in the ModelList
to compute Gₛ.
We can use the stomatal conductance model of Medlyn et al. (2011) as an example to compute it. It is implemented with the Medlyn
structure. We can then create a ModelList
with the two models:
ModelList(Fvcb(), Medlyn(0.03, 12.0))
PlantSimEngine.DependencyGraph{Dict{Symbol, PlantSimEngine.SoftDependencyNode}}(Dict{Symbol, PlantSimEngine.SoftDependencyNode}(:photosynthesis => Fvcb{Float64}
), Dict{Symbol, DataType}())TimeStepTable{Status{(:aPPFD, :Tₗ, :Cₛ,...}(1 x 7):
╭─────┬─────────┬─────────┬─────────┬─────────┬─────────┬─────────┬─────────╮
│ Row │ aPPFD │ Tₗ │ Cₛ │ A │ Gₛ │ Cᵢ │ Dₗ │
│ │ Float64 │ Float64 │ Float64 │ Float64 │ Float64 │ Float64 │ Float64 │
├─────┼─────────┼─────────┼─────────┼─────────┼─────────┼─────────┼─────────┤
│ 1 │ -Inf │ -Inf │ -Inf │ -Inf │ -Inf │ -Inf │ -Inf │
╰─────┴─────────┴─────────┴─────────┴─────────┴─────────┴─────────┴─────────╯
Now this instantiation returns some warnings saying we need to initialize some variables.
The Fvcb
model requires the following variables as inputs:
inputs(Fvcb())
(:aPPFD, :Tₗ, :Cₛ)
And the Medlyn
model requires the following variables:
inputs(Medlyn(0.03, 12.0))
(:Dₗ, :Cₛ, :A)
We see that A
is needed as input of Medlyn
, but we also know that it is an output of Fvcb
. This is why we prefer using to_initialize
from PlantSimEngine.jl
instead of inputs
, because it returns only the variables that need to be initialized, considering that some inputs are duplicated between models, and some are computed by other models (they are outputs of a model):
to_initialize(ModelList(Fvcb(), Medlyn(0.03, 12.0)))
(photosynthesis = (:aPPFD, :Tₗ, :Cₛ), stomatal_conductance = (:Dₗ, :Cₛ))
The most straightforward way of initializing a model list is by giving the initializations to the status
keyword argument during instantiation:
m = ModelList(
Fvcb(),
Medlyn(0.03, 12.0),
status = (Tₗ = 25.0, aPPFD = 1000.0, Cₛ = 400.0, Dₗ = 0.82)
)
PlantSimEngine.DependencyGraph{Dict{Symbol, PlantSimEngine.SoftDependencyNode}}(Dict{Symbol, PlantSimEngine.SoftDependencyNode}(:photosynthesis => Fvcb{Float64}
), Dict{Symbol, DataType}())TimeStepTable{Status{(:aPPFD, :Tₗ, :Cₛ,...}(1 x 7):
╭─────┬─────────┬─────────┬─────────┬─────────┬─────────┬─────────┬─────────╮
│ Row │ aPPFD │ Tₗ │ Cₛ │ A │ Gₛ │ Cᵢ │ Dₗ │
│ │ Float64 │ Float64 │ Float64 │ Float64 │ Float64 │ Float64 │ Float64 │
├─────┼─────────┼─────────┼─────────┼─────────┼─────────┼─────────┼─────────┤
│ 1 │ 1000.0 │ 25.0 │ 400.0 │ -Inf │ -Inf │ -Inf │ 0.82 │
╰─────┴─────────┴─────────┴─────────┴─────────┴─────────┴─────────┴─────────╯
Our component models structure is now fully parameterized and initialized for a simulation!
Let's simulate it:
run!(m)
The models included in the package are listed in their own section, i.e. here for photosynthesis.