# PlantBiophysics.jl: a high-performance, modular software for prototyping and scaling biophysical models from leaf to canopy

Rémi Vezy, Simon Treillou, Samuel Mackeown, Emilie Peynaud, Raphaël P.A. Perez, Thomas Arsouze, Jean Dauzat

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Rémi Vezy, Simon Treillou, Samuel Mackeown, Emilie Peynaud, Raphaël P.A. Perez, et al., Plant-Biophysics.jl: a high-performance, modular software for prototyping and scaling biophysical models from leaf to canopy. in *silico* Plants, 2025, 7 (2), pp.diaf021. 10.1093/insilicoplants/diaf021. *hal-05376887*

HAL Id: hal-05376887

https://hal.inrae.fr/hal-05376887v1

Submitted on 21 Nov 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
# Methods for Plant Modeling

## PlantBiophysics.jl: a high-performance, modular software for prototyping and scaling biophysical models from leaf to canopy

Rémi Vezy<fnref n="1" />,<fnref n="2" />,*† ID, Simon Treillou<fnref n="1" />,<fnref n="2" />,† ID, Samuel Mackeown<fnref n="1" />,<fnref n="2" />, Emilie Peynaud<fnref n="1" />,<fnref n="2" /> ID,
Raphaël P.A. Perez<fnref n="3" /> ID, Thomas Arsouze<fnref n="1" />,<fnref n="2" /> ID, Jean Dauzat<fnref n="1" />,<fnref n="2" />

Handling Editor: Steve Long

**Abstract.** Process-based biophysical models are essential for understanding and predicting plant responses to environment, yet their development and application remains hindered by the 'two-language problem': models are typically prototyped in interactive, slow languages and then reimplemented in compiled, fast languages for simulation. This workflow requires more human resources, and it may introduce delays, errors, and complexity for model evaluation and prototyping iteration. We present PlantBiophysics.jl, a pure Julia package that unifies model prototyping, calibration, simulation, and evaluation of biophysical models within a single, high-performance and interactive environment. The package currently supports four key processes: light interception, energy balance, photosynthesis, and stomatal conductance. It provides automated model composition and coupling through a modular and extensible architecture. The package supports uncertainty propagation and compatibility with 3D plant structures and dynamic meteorological inputs. We benchmarked PlantBiophysics.jl against the widely used R package plantecophys and the Julia package LeafGasExchange.jl. Simulations ran up to ~6700× and ~460× faster, respectively, while improving predictive accuracy (e.g. net assimilation nRMSE: 5% vs. 9% and 12%, respectively). By allowing models to be easily swapped, coupled, or extended, and by supporting simulations at scales ranging from single leaves to 3D plant canopies, PlantBiophysics.jl provides a robust, efficient, and transparent tool for plant ecophysiologists and modellers.

**Keywords:** energy balance; photosynthesis; stomatal conductance; gas exchange; biophysics; FSPM

## 1. INTRODUCTION

Understanding and predicting how plants exchange energy, water, and carbon with their environment requires accurate modelling of biophysical processes, including photosynthesis, stomatal conductance, energy balance, and evapotranspiration. These processes span scales from the individual organ to the whole canopy and are tightly influenced by environmental drivers such as light, air temperature, humidity, and wind, creating a highly dynamic and interconnected system (Monteith and Unsworth 2013, Gutschick 2016). Because they form the foundation of the soil-plant-atmosphere continuum, their accurate representation is crucial across a wide range of disciplines, including plant science, agronomy, ecology and climate research (Evers et al. 2018, Wang et al. 2021).

Biophysical models serve as a bridge between detailed physiological knowledge and large-scale predictions. They enable researchers to extrapolate local leaf-level measurements to canopy or ecosystem scale, assess plant responses to environmental or management changes, optimize genotype selection and practice combinations, and study emergent properties in complex plant-environment systems. Early applications, such as the work by Parkhurst and Loucks (1972), already demonstrated their value by showing how simple leaf energy balance models could explain the adaptive significance of leaf size under different climates. Today, modules for photosynthesis and energy balance are routinely embedded in global Earth system models (Lombardozzi et al. 2015, Rogers et al. 2017, Riley et al. 2018), and continue to be central to studies of climate

<footnote n="1">¹CIRAD, UMR AMAP, Montpellier F-34398, France</footnote>

<footnote n="2">²AMAP, Univ Montpellier, CIRAD, CNRS, INRAE, IRD, Montpellier, France</footnote>

<footnote n="3">³CIRAD, UMR AGAP Institute, Montpellier F-34398, France</footnote>

<footnote>*Corresponding author. UMR AMAP, CIRAD - TA A51/PS2 34398 Montpellier cedex 5, France. E-mail: remi.vezy@cirad.fr</footnote>

<footnote>†These authors should be considered joint first author.</footnote>
change impacts on vegetation function and productivity (Crous 2019, Vezy et al. 2020).

Yet despite decades of research, several limitations still constrain the development and use of biophysical models. One challenge is the proliferation of alternative formulations for the same process. Stomatal conductance, for example, can be represented by multiple models that differ in their scale, assumptions, formulation, or empirical corrections (Damour et al. 2010, Buckley 2017). Even minor corrections, such as the correction proposed by Schymanski and Or (2017) to the Penman-Monteith model (Monteith and Unsworth 2013) can significantly improve accuracy while requiring minimal alterations in the mathematical formulation. These alternative implementations reflect both progress in scientific understanding and differences in modelling objectives, whether related to simulation scale, data availability, or computational demands.

A second challenge is that combining models into a single coherent workflow remains non-trivial. Different formalisms may rely on incompatible inputs or require specific model call sequences, making integration of new models laborious. Many models are written in compiled languages such as C++, Fortran, or Java, and are embedded in complex architectures where modifying or swapping individual components is difficult. Adding new models or adjusting the model chain often requires a deep understanding of the underlying software and considerable effort for revalidation.

A third and persistent issue is the trade-off between computational performance and usability. High-performance compiled models like STICS (Brisson et al. 1998) or MAESPA (Duursma and Medlyn 2012) are suitable for large-scale simulations but tend to be monolithic and inflexible. By contrast, packages written in interpreted languages such as R are easier to use and modify, and have gained popularity for prototyping. For example, the plantecophys package (Duursma 2015) is widely used in the ecophysiology community, alongside other tools like tealeaves (Muir 2019) and photosynthesis (Stinziano et al. 2021). However, these tools often lack the computational speed needed for large-scale or high-frequency simulations (Liu et al. 2021), particularly when modelling complex canopies with millions of leaves or long time series.

This workflow gap, often referred to as the 'two-language problem', means that models must be written twice: once in a flexible, slow language for development and prototyping, and again in a fast, compiled language for deployment. This practice introduces redundancy, slows iteration, and increases the likelihood of implementation errors.

To bridge this gap, plant modellers increasingly need tools that combine speed, modularity, and ease of use. An ideal platform would support rapid prototyping, flexible coupling of components, and large-scale simulation—all within a single, high-performance environment. It would also provide transparent, well-tested implementations, and facilitate tasks like uncertainty analysis, model comparison, and calibration.

Here, we introduce PlantBiophysics.jl, a pure Julia package implemented using the PlantSimEngine.jl modelling framework (Vezy 2023). Julia is a modern programming language that combines the speed of compiled languages with the ease and flexibility of scripting languages (Bezanson et al. 2017). It is increasingly used in biological modelling (Roesch et al.

2023), notably within the functional-structural plant modelling community through the Virtual Plant Laboratory (Morales et al. 2025) and its associated ecosystem of packages (VPLverse), which includes PlantBiophysics.jl. Built on a modular architecture, PlantBiophysics.jl supports biophysical model prototyping, calibration, and simulation in a single high-performance environment. It currently includes four key processes: light interception, energy balance, photosynthesis, and stomatal conductance, with automatic model coupling and support for uncertainty propagation. Model coupling is handled automatically by PlantSimEngine.jl, and because the framework allows passing arbitrary numerical types, the package also supports advanced use cases such as uncertainty propagation. Simulations can be performed on a single leaf, across many leaves, or using a multiscale structural representation. The package is available from the general Julia registry under the MIT open-source license and follows community standards best practices (Stinziano et al. 2021). By bridging the gap between ease of use and computational efficiency, it enables simulations that scale from individual leaves to whole-plant 3D scenes without requiring model re-implementation.

In this article we describe the core design of PlantBiophysics.jl, explain how its modular system enables efficient model coupling, present an evaluation on data, and compare its accuracy and computational performance with two widely used packages, plantecophys (Duursma 2015) and LeafGasExchange.jl (Yun et al. 2020). We demonstrate that PlantBiophysics.jl provides a reliable, modular, accurate, and high-performance platform that lowers the barrier to developing, fitting and simulating coupled biophysical processes across scales.

## 2. MODEL AND DATA DESCRIPTION

### 2.1 Package design, architecture and modularity

PlantBiophysics.jl is a Julia package designed to prototype, calibrate and simulate biophysical models efficiently and reproducibly. It is built on the PlantSimEngine.jl framework, with optional support from PlantMeteo.jl and MultiScaleTreeGraph.jl. PlantSimEngine.jl (Vezy 2023) provides standardized tools for model development, automatically handling simulation workflows and model coupling. PlantMeteo.jl manages meteorological inputs and is particularly useful when higher performance is required compared with a standard 'DataFrame'. MultiScale TreeGraph.jl (Vezy 2025) offers data structures and algorithms for an optional multiscale representation of plant scenes, from individual organs to entire canopies.

This architecture separates structural description of the system (scene, plants, organs) from its functional processes (e.g. photosynthesis, conductance, energy balance). It ensures that each model is implemented as an independent component, can be unit-tested in isolation, and can be flexibly replaced, extended or coupled. This design not only facilitates model development but also supports model comparison and coupling.

PlantBiophysics.jl can operate at any scale, from single-leaf calculations to simulations of entire 3D scenes. For simple applications such as model calibration, users can run the package on a single object (e.g. a leaf with gas-exchange measurements). For more complex simulations, the scene can contain arbitrary
objects such as plants, soil layers, or artificial structures. Each
object can be decomposed into components at multiple scales.
For example, a tree can be described as a single ‘big leaf’, as a
collection of axes, or as a fully resolved set of internodes and
leaves (Godin and Caraglio 1998). Then, a model is associated
to each scale, either simulating the process or integrating the
results at a coarser scale. Further details are provided in
MultiScaleTreeGraph.jl documentation (Vezy 2025).

During a simulation, the functional state of each component
evolves according to the selected models and the current envir-
onmental conditions. Measured variables (e.g. air temperature,
radiation) can be forced as inputs, while other state variables
are computed by models. PlantSimEngine.jl automatically
builds a dependency graph between models, ensuring they
are executed in the correct order and allowing consistent coup-
ling as well as parallel computation across time steps and ob-
jects (Vezy 2023).

At the implementation level, PlantBiophysics.jl leverages
Julia's multiple dispatch and PlantSimEngine.jl's type-stable
architecture to define modular model components that interact
through a minimal interface. Each model needs to implement
only three methods (i.e. 'inputs', 'outputs', and 'run!') to be
compatible with the simulation engine. This approach makes
it straightforward to swap model formulations without altering
the rest of the simulation, e.g. replacing a vapour-pressure def-
icit (VPD)-based stomatal conductance model (Medlyn et al.
2011) with a water potential-based alternative (Tuzet et al.
2003). The execution order and data dependencies are auto-
matically managed by the dependency graph system in
PlantSimEngine.jl, which reduces the risk of integration errors.

A further advantage of this modular structure is type-generic
design. All code is written generically with respect to numerical
types, enabling the use of any standard types (e.g. 'Float64'),
but also units or uncertainty-aware types such as 'Particles'
from MonteCarloMeasurements.jl (Carlson 2020) or
Measurements.jl (Giordano 2016). This flexibility allows

advanced features such as uncertainty propagation or sensitivity
analysis to be introduced simply by changing the input type,
without modifying model equations or coupling logic.

This design philosophy ensures that new models can be im-
plemented with minimal effort, and that performance improve-
ments to the framework automatically benefit all user-defined
models. By clearly separating model logic from computational
infrastructure, the framework reduces long-term maintenance
overhead and allows users to focus on model development.
New biophysical processes can be added as independent com-
ponents and seamlessly combined with existing ones, making
PlantBiophysics.jl particularly well suited for model intercom-
parison, exploratory simulations, and multiscale modelling
workflows.

## 2.2 Model instructions

2.2.1 *Model implementation.* Implementing a new model in PlantBiophysics.jl is straightforward thanks to the PlantSimEngine.jl framework. PlantSimEngine.jl provides high-level methods that automatically handle model coupling via a dependency graph, along with built-in checks, memory pre-allocation, and parallelization over time steps and objects to improve efficiency.

As an example, the Beer–Lambert law of light extinction can be
implemented in 14 lines of code (Fig. 1). First, a structure is
declared to name the model and list its parameter values. This
structure must be a subtype of the abstract type corresponding
to its process, which in this case is ‘AbstractLight_
InterceptionModel’ (a full list of available abstract types is pro-
vided in the documentation: https://vezy.github.io/Plant
Biophysics.jl/stable/extending/implement_a_model). Using a
parametric type (i.e. the {T} after the structure name) is recom-
mended to ensure type stability and composability (see https://
docs.julialang.org/en/v1/manual/types/#Parametric-Types). This
approach allows the model to accept any data type, such as

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Code snippet showing an implementation of the Beer–Lambert law of light extinction in PlantBiophysics.jl using PlantSimEngine.jl. Please note that the syntax is subject to changes with new versions of PlantSimEngine.jl. Always refer to the online documentation of the stable release here: https://virtualplantlab.github.io/PlantSimEngine.jl.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Code snippet of a model coupling into a 'ModelList': using the 'Monteith' model for the energy balance process (Monteith and Unsworth 2013), 'Fvcb' for photosynthesis (Farquhar et al. 1980), and 'Medlyn' for stomatal conductance (Medlyn et al. 2011). The 'ModelList' takes a list of models and automatically couples the models using a directed acyclic dependency graph. It also takes a named argument called 'status' to initialize variables that are not computed by any model. The list of variables requiring initialization for a given 'ModelList' can be obtained by applying the 'to_initialize' function to the 'ModelList'. 'Ra_SW_f': absorbed shortwave radiation flux (W m⁻²), 'sky: sky view factor for both faces of the object (0–2), 'aPPFD': absorbed photosynthetic photon flux density (µmol photon m⁻² s⁻¹), 'd': characteristic dimension (e.g. leaf width, m).</figcaption></figure>

units (see e.g. the Unitful.jl package) or 'Particles' for error propagation, while preserving the type stability needed for high performance.

The model's algorithm is then implemented in a new method for the 'run!' function from PlantSimEngine.jl. Following Julia convention, the exclamation mark indicates that the function modifies its inputs. Finally, two methods for the 'inputs_' and 'outputs_' functions specify the model's required inputs and outputs. These functions, which follow Julia's convention of using an underscore suffix for internal methods, ensure that models are correctly initialized within a simulation. The function names are also prefixed with 'PlantSimEngine'. to make it clear they extend existing functions from the framework rather than creating new ones.

With these few steps, the Beer–Lambert model can be used either in isolation or coupled with other models. This flexibility is especially valuable for non-expert users, who can integrate their models into a broader simulation framework alongside pre-existing components.

### 2.2.2 Model coupling.
Model coupling in PlantBiophysics.jl is managed automatically via the 'ModelList' structure from PlantSimEngine.jl. To couple multiple biophysical processes, users simply list the desired models within a 'ModelList' (Fig. 2). PlantSimEngine.jl then determines the correct execution order using a (multi-scale) directed acyclic dependency graph, which resolves data flow between model inputs and outputs (soft-coupling) and supports hard-coupling when a model explicitly calls another, such as during iterative computations.

This design allows users to focus on 'what' processes to include and 'which' models to use rather than 'how' to couple them. The 'ModelList' can also be initialized with a 'status' argument, providing initial values for variables that are not computed by any model. The list of variables requiring user initialization can be obtained with the 'to_initialize' function. For multi-scale simulations on plant graphs (e.g. using MultiScaleTreeGraph.jl),

models can be assigned to specific spatial scales through a mapping dictionary. Although this more advanced configuration is not shown here, detailed examples are available in the PlantBiophysics.jl and PlantSimEngine.jl documentation.

### 2.2.3 Alternative model formulations.
Swapping between alternative model formulations is equally straightforward. To test the effect of a different process representation, the user simply replaces the model's name in the 'ModelList' while keeping the rest of the setup unchanged. For instance, 'Medlyn' can be replaced with 'Tuzet' to compute stomatal conductance to CO₂ using Tuzet's model, while still coupling it with the same photosynthesis and energy balance models (Fig. 3).

This capability supports rapid comparison of alternative process representations in isolation or in coupled simulations. It can also be used to replace a process model with one that forces a variable to match observations, thereby reducing the degrees of freedom in a larger model. Because PlantSimEngine.jl automatically updates the dependency graph, these changes require no additional adjustments to the simulation code.

### 2.2.4 Model fitting and file IO.
PlantBiophysics.jl provides dedicated methods for parameter estimation through the generic 'fit' function from PlantSimEngine.jl. This functionality allows users to calibrate a model directly from empirical data, facilitating integration of field or laboratory measurements into the modelling workflow.

The specific arguments accepted by 'fit' depend on the model and the fitting method, but typically include the model type (e.g. 'Fvcb'), a 'DataFrame' containing the experimental data and optional keyword arguments, such as fixed parameter values or environmental conditions to further guide the fitting procedure (see Fig. 4 for an example using 'Fvcb').

A common challenge in model calibration is that gas-exchange data are stored in different, often inconsistent formats. Each instrument vendor (e.g. LI-COR, WALZ) uses its
<figure><img src="image_5.png" /></figure>

```julia
models = ModelList(
    Monteith(),
    Fvcb(VcMaxRef=200.0, JMaxRef= 250.0),
    Tuzet(g0=0.03, g1=12.0, Ψv=-1.5, sf=2.0, Γ=30.0),
    status=(Ra_SW_f=13.747, skyfraction=1.0, aPPFD=1500.0, d=0.03, Ψl=-1.0)
)
```

**Figure 3.** Code snippet of a model coupling as in Fig. 2, but using the model of Tuzet et al. (2003) instead of the model of Medlyn et al. (2015) to compute stomatal conductance to CO₂. 'Ra_SW_f': absorbed shortwave radiation flux (W m⁻²), 'skyfraction': view factor between the object and the sky for both faces (0–2), 'aPPFD': absorbed photosynthetic photon flux density (µmol photon m⁻² s⁻¹), 'd': characteristic dimension (e.g. leaf width, m), 'Ψl': leaf water potential (MPa).

<figure><img src="image_6.png" /></figure>

```julia
df = read_licor6800("file.csv")

VcMaxRef, JMaxRef, RdRef, TPURef = fit(Fvcb, df; T_r=25.0)
```

**Figure 4.** Code snippet showing how to import A–Cᵢ response curve data measured with a LI-COR 6800 portable gas-exchange system and calibrate the Fvcb model (Farquhar et al. 1980, von Caemmerer and Farquhar 1981).

own file structure and variable naming conventions, and even within a given vendor's format, users may modify the format to fit their needs. As a result, there is no single 'standard' way to obtain ready-to-use data from the instrument.

To address this, PlantBiophysics.jl provides dedicated functions for reading data from the most common devices: 'read_licor6400' and 'read_licor6800' for LI-COR instruments, 'read_walz' for the WALZ GFS-3000, and 'read_ciras4' for PPSystem CIRAS-4. These functions automatically parse the vendor-specific formats and return a standardized 'Frame' with columns such as Tᵢ (leaf temperature, °C), aPPFD (absorbed photosynthetic photon flux density, µmol m⁻² s⁻¹), Cᵢ (intercellular CO₂ concentration, µmol mol⁻¹), A (net assimilation rate, µmol m⁻² s⁻¹), Cₐ (ambient CO₂ concentration, µmol mol⁻¹), Dᵢ (leaf-to-air vapour pressure deficit, kPa), and Gₛ (stomatal conductance to water vapour, mol m⁻² s⁻¹), ready for use in the fitting routines.

For users seeking a reproducible, community-agreed format, PlantBiophysics.jl also supports the ESS-DIVE leaf-level gas-exchange standard (Ely et al. 2021) via the 'read_ess_dive' function. This requires that users first prepare their data to match the ESS-DIVE specification, after which the function can directly ingest the file into a compatible 'Frame'.

### 2.2.5 Step-by-step example of a simulation.
Running a simulation in PlantBiophysics.jl requires only two main elements: a list of models describing the biophysical processes to simulate (the 'ModelList'), and a set of meteorological conditions. Simulations are launched using the 'run!' function, which automatically handles model execution based on the dependencies defines the model list.

In its simplest form, a simulation can be performed on a single leaf for a single time step (Fig. 5). The user defines a 'ModelList' containing the desired models, such as those for energy balance, photosynthesis, and stomatal conductance, and provides a single set of meteorological conditions as an 'Atmosphere' provided by PlantMeteo.jl, which pre-computes atmosphere-related variables based on user inputs (i.e. vapour pressure, vapour pressure deficit, etc...). The run! function then computes all outputs associated with the selected processes, updating the system's state accordingly.

For simulations over multiple time steps, the meteorological input can be provided in any format that implements the Tables.jl interface. This includes common data structures such as DataFrames, vectors of named tuples, or even more complex sources such as SQLite databases (see e.g. SQLite.jl) or TimeStepTables (from PlantMeteo.jl). The
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Code snippet showing an example simulation of coupled energy balance, photosynthesis, and stomatal conductance processes for a single leaf under a specified meteorological condition.</figcaption></figure>

'run!' function automatically iterates over all time steps, ensuring that the models are executed in the right order and consistently over the period.

Simulating multiple leaves or other objects follows the same principles. A vector of 'ModelList's can be passed to the 'run!' function, each corresponding to the parameterization of a different object. Alternatively, for more complex structures such as plants or scenes, a MultiScaleTreeGraph can be provided. This allows simulations across multiple spatial scales, with each node (e.g. leaf, branch, or entire plant) associated with its own set of models. In both cases, the simulation engine ensures that each model is evaluated at the proper time depending on its input and output variables thanks to the (multiscale-) dependency graph.

## 2.3 Uncertainty propagation

In most ecophysiological and biophysical modelling frameworks, input data such as meteorological variables, morphological traits, or model parameters are treated as fixed scalar values. However, these variables are almost always subject to measurement errors, instrument noise, or biological variability across individuals or conditions. While these uncertainties are often acknowledged qualitatively, they are rarely formally propagated through the model equations. As a result, model outputs are typically reported without any assessment of the reliability or robustness of the predictions.

PlantBiophysics.jl addresses this gap by natively supporting uncertainty propagation as a first-class feature of the modelling framework (like any other feature implemented via types, e.g. units). Thanks to its foundation in Julia and the design principles of PlantSimEngine.jl, all models in the package are fully type-stable and composable, allowing them to operate not only on scalars but also on distributions or uncertainty-aware numeric types.

In particular, PlantBiophysics.jl integrates seamlessly with Measurements.jl (Giordano 2016) and MonteCarlo Measurements.jl (Carlson 2020), which provide tools for linear and Monte Carlo-based uncertainty propagation. Under this approach, scalar inputs are replaced with a special type: a Measurement or a Particle, respectively. The Measurement type defines a nominal value and a standard deviation, that are defined

independently on the inputs and propagated by operator overloading. The Particle type represents arbitrary distributions as a vector of unweighted particles, enabling accurate representation and efficient propagation of nonlinear uncertainty. For example, a normal distribution for radiation with 10% standard deviation, or a uniform distribution for leaf characteristic length. Users of PlantBiophysics.jl are encouraged to use the Monte Carlo-based uncertainty propagation because most algorithms in ecophysiology are highly non-linear. In this set-up, simulations are performed independently across random samples (i.e. particles), automatically yielding output distributions that reflect how input uncertainty affects model predictions. The method is fully compatible with nonlinear models and provides a practical way to quantify confidence intervals around simulated variables. The accuracy of the propagated uncertainty depends on the number of samples (n), with higher values offering smoother and more reliable distributions at the cost of increased computation time.

Uncertainty propagation can be applied to any part of the simulation, from energy balance and leaf temperature to photosynthesis and stomatal conductance. It is especially useful for analysing measurement-derived inputs such as radiation, wind speed, or anatomical traits. The resulting outputs can be evaluated using classical performance metrics (e.g. modelling efficiency or bias), with the added benefit of accessing the full distributions instead of point estimates.

While powerful, Monte Carlo propagation requires some care. For instance, some variables such as characteristic length or radiation must remain strictly positive, and inappropriate choices of distributions (e.g. symmetric Gaussians centred near zero) can lead to invalid or non-physical simulations. Users are encouraged to select uncertainty distributions that reflect the physical nature and empirical variability of each parameter.

## 2.4 Implemented processes and models

PlantBiophysics.jl currently implements four core biophysical processes: light interception, energy balance, photosynthesis, and stomatal conductance (Table 1). Thanks to PlantSimEngine.jl, each process can be simulated with one or more alternative model formulations. These models are independent, composable, and can be easily substituted or extended
<table><caption>Table 1. List of models for the four different processes in PlantBiophysics.jl available at the time of writing.</caption><thead><tr><th>Process</th><th>Function name</th><th>Model</th><th>Description</th></tr></thead><tbody><tr><td>Light interception</td><td>light_interception</td><td>Beer</td><td>The Beer-Lambert law of light extinction</td></tr><tr><td>Energy balance (latent and sensible fluxes)</td><td>energy_balance</td><td>Monteith</td><td>Monteith and Unsworth (2013) model</td></tr><tr><td rowspan="5">Photosynthesis</td><td rowspan="5">Photosynthesis</td><td>FvcbRaw</td><td>Farquhar-von Caemmerer-Berry (FvCB) model for C3 photosynthesis (Farquhar et al. 1980, von Caemmerer and Farquhar 1981)</td></tr><tr><td>Fvcb</td><td>FvcbRaw coupled with a stomatal conductance model with an analytical resolution (Baldocchi 1994).</td></tr><tr><td>FvcbIter</td><td>Fvcb with an iterative resolution</td></tr><tr><td>ConstantA</td><td>Constant assimilation</td></tr><tr><td>ConstantAGs</td><td>Constant assimilation, coupled to a stomatal conductance model to compute C<sub>i</sub></td></tr><tr><td rowspan="3">Stomatal conductance</td><td rowspan="3">stomatal_conductance</td><td>Medlyn</td><td>Medlyn et al. (2011) model</td></tr><tr><td>Tuzet</td><td>Tuzet et al. (2003) model</td></tr><tr><td>ConstantGs</td><td>Constant stomatal conductance</td></tr></tbody></table>

without changing the overall simulation workflow. This modularity facilitates both model intercomparison and accurate reproduction of classical formalisms.

### 2.4.1 Light interception.
Light interception is computed using the Beer-Lambert law of light extinction. This process calculates the fraction of incident radiation absorbed by an object based on the leaf area index and is often used for homogeneous canopies (Ponce de León and Bailey 2019, Bailey et al. 2020).

### 2.4.2 Photosynthesis.
Photosynthesis is simulated according to the Farquhar-von Caemmerer-Berry C3 photosynthesis model (Farquhar et al. 1980, von Caemmerer and Farquhar 1981), also known as FvCB model, defining the rate of photosynthesis $A_n$ as the smallest of three potential rates of carboxylation minus a dark respiration rate term $R_d$. The model main equation can be written as:

$$A_n = \min(W_j, W_v, W_p) - R_d \quad (1)$$

with $A_n$ the net carbon assimilation; $W_j$ the electron-transport-limited rate of CO₂ assimilation; $W_v$ the Rubisco-carboxylation-limited rate, i.e. when the kinetics of the RuBisCo enzyme for CO₂ is at its maximum capacity, which happens mostly when the CO₂ concentration in the stomata is too low; $W_p$ the Triose Phosphate Utilization (TPU)-limited rate, and $R_d$ the rate of mitochondrial respiration at the given leaf temperature. All these variables are expressed in µmol CO₂ m⁻² s⁻¹.

The electron-transport-limited rate $W_j$ is computed using $V_j$ as follows:

$$W_j = V_j \left( \frac{C_i - \Gamma^s}{C_i + 2\Gamma^s} \right) \quad (2)$$

with $V_j = \frac{1}{4\cdot2\theta}\alpha \cdot aPPFD + J_{Max} - \sqrt{(\alpha \cdot aPPFD + J_{Max})^2 - 4\alpha \cdot \theta \cdot J_{Max} \cdot aPPFD}$),
$C_i$ the intercellular CO₂ concentration (ppm), $\Gamma^s$ the CO₂ compensation point (µmol mol⁻¹), $\theta$ the parameter determining the shape of the non-rectangular hyperbola, $\alpha$ the quantum yield of electron transport (mol_e mol⁻¹ quanta), aPPFD the absorbed photosynthetic photon flux density

($\mu\text{mol quanta m}^{-2}\text{s}^{-1}$) and $J_{Max}$ the maximum rate of electron transport derived from the reference term $J_{MaxRef}$ using the Arrhenius function ($\mu\text{mol CO}_2\text{m}^{-2}\text{s}^{-1}$). The Rubisco-carboxylation-limited rate is given by:

$$W_v = V_{cMax} \left( \frac{C_i - \Gamma^s}{C_i + K_m} \right) \quad (3)$$

with $V_{cMax}$ derived from the reference term $V_{cMaxRef}$ using the Arrhenius function ($\mu\text{mol CO}_2\text{m}^{-2}\text{s}^{-1}$) and $K_m$ the effective Michaelis-Menten coefficient for CO₂. And the TPU-limited rate is:

$$W_p = \frac{3 \cdot \text{TPU}(C_i - \Gamma^s)}{C_i - (1 + 3a_g) \cdot \Gamma^s} \quad (4)$$

simplified in $W_p = 3 \times \text{TPU}$ when $a_g$ is 0 (Lombardozzi et al. 2018), $a_g$ being a unitless scalar between 0 and 1 based on glycolate recycling.

PlantBiophysics.jl has three different implementations of photosynthesis. 'FvCBRaw' is the implementation of the FvCB model as is, i.e. without any coupling with a conductance model. All variables needed for the 'FvCBRaw' model are direct inputs. 'FvCBIter' must be used in conjunction with a stomatal conductance model. Given a CO₂ boundary conductance $G_{bc}$, the algorithm iterates until a stable solution for the couple ($G_s, A_n$) is found. Finally, 'FvCB' is an analytical resolution of the stomatal conductance-photosynthesis coupling following Baldocchi (1994) and requiring the stomatal atmospheric CO₂ concentration $C_s$ as an input.

### 2.4.3 Stomatal conductance.
Three models have been implemented to model stomatal conductance to CO₂: a constant stomatal conductance model, the model presented in Medlyn et al. (2011), and the model presented in Tuzet et al. (2003), herein called the Medlyn and Tuzet models, respectively. Medlyn's model is as follows:

$$G_s = g_0 + \left( 1 + \frac{g_1}{\sqrt{D_l}} \right) \frac{A_n}{C_s} \quad (5)$$

where $g_0$ and $g_1$ are two parameters, $D_l$ the water vapour pressure deficit between the leaf and the ambient air (kPa), $A_n$ the
net assimilation of carbon dioxide in µmolCO₂m⁻²s⁻¹, and C₅ the CO₂ concentration at the leaf surface (ppm).

The original Medlyn model uses the air CO₂ concentration (Cₐ) rather than the leaf-surface concentration (C₅), likely reflecting its derivation from portable gas-exchange measurements where air is strongly mixed and the boundary layer around the leaf can be considered negligibly thin, making Cₐ ≈ C₅. In natural conditions, this assumption does not always hold, as boundary layers may create a measurable difference between the two concentrations. PlantBiophysics.jl therefore uses C₅ instead of Cₐ to provide a more general formulation. Users can either supply C₅ from a boundary-layer model or simply consider that C₅ = Cₐ for matching with the original model description.

The model of Tuzet is defined as:

$$G_s = g_0 + g_1 \cdot \frac{A_n}{C_s - \Gamma} \cdot f_{\Psi_v} \quad (6)$$

where g₀ and g₁ are two parameters, Γ is the CO₂ compensation point (mol mol⁻¹), and fΨ_v is defined as:

$$f_{\Psi_v} = \frac{1 + e^{s_f \Psi_f}}{1 + e^{s_f (\Psi_f - \Psi_v)}} \quad (7)$$

where Ψ_v the bulk leaf water potential, Ψ_f is a reference water potential and s_f is a stomatal closure sensitivity to water potential.

While the original Tuzet et al. (2003) model used Cᵢ, PlantBiophysics.jl uses C₅ for consistency with the other implementations and because model parameters are typically fitted from gas-exchange data where C₅ is the directly available quantity. This choice also avoids introducing extra coupling complexity between Aₙ, G₅ and Cᵢ that are interdependent.

### 2.4.4 Energy balance.
Leaves play a critical role in the exchange of energy between plants and their environment, and a key challenge is to achieve an accurate balance between the energy gained, lost, and stored by all leaves at any time. The classical energy balance equation described by Monteith and Unsworth (2013) can be used to calculate the energy balance of an object:

$$R_n + M = H + \lambda E + G \quad (8)$$

where Rₙ is the energy flux density absorbed by the object via radiation, M is the metabolic energy flux density, λE is the latent heat flux density (i.e. the energy needed to change the state from e.g. vapour to liquid water), H is the sensible heat flux density, and G the conductively transmitted heat flux density. All quantities are usually expressed in W m⁻².

In the case of a leaf, the energy balance equation can be further simplified by making several approximations:

* The leaf is assumed to be infinitely thin, which implies that the site of latent energy exchange (the mesophyll) and the site of sensible energy exchange (the epidermis and cuticle) are at the same average temperature. As a result, the leaf can be treated as a single temperature object with no temperature gradient.

* Over a sufficiently long period, there is no significant net energy storage within the leaf. This allows us to neglect G in the equation.

* The metabolic energy (M) contribution to the overall energy balance is considered negligible compared with other terms, and can be neglected in the equation.

Under these hypotheses, PlantBiophysics.jl uses a simplified energy balance written as:

$$R_n \approx H + \lambda E \quad (9)$$

### 2.4.5 Coupled photosynthesis, stomatal conductance, and energy balance.
PlantBiophysics.jl implements biophysics models using the PlantSimEngine.jl simulation framework, together providing an integrated system for simulating the tight coupling between photosynthesis, stomatal conductance, and leaf energy balance. These three processes are dynamically interdependent: leaf temperature influences photosynthesis and transpiration; transpiration affects stomatal conductance; and conductance regulates CO₂ diffusion, closing the feedback loop. This coupling mechanism is illustrated in Fig. 6 with the example of using the FvCB, Medlyn and Monteith models.

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Diagram of an example model coupling in PlantBiophysics.jl used to compute energy balance, photosynthesis, and stomatal conductance for a C3 plant. Photosynthesis and stomatal conductance are coupled using an analytical resolution, and energy balance is computed iteratively based on a convergence criterion over leaf temperature.</figcaption></figure>
**2.4.5.1 Coupling photosynthesis and stomatal conductance.** Net CO₂ assimilation ($A_n$) depends on the intercellular CO₂ concentration ($C_i$), which is regulated by stomatal conductance ($G_s$). Conversely, in Medlyn's and Tuzet's models, $G_s$ is partially controlled by $A_n$, forming a bidirectional dependency. In PlantBiophysics.jl, this feedback is resolved either analytically using the Baldocchi (1994) formulation (available via the 'FvCB' model), or iteratively using 'FvCBIter', which iterates until a stable $A_n-C_i$ solution is found.

**2.4.5.2 Coupling gas exchange with energy balance.** The energy balance determines leaf temperature ($T_l$), which directly affects photosynthesis rate through enzyme kinetics and electron transport, and stomatal conductance by modulating the leaf-to-air vapour pressure deficit ($D_l$). In turn, stomatal conductance influences transpiration and thus the latent heat flux ($\lambda E$), affecting leaf temperature. This creates a second feedback loop between gas exchange and thermodynamics.

PlantBiophysics.jl resolves this interaction using a fixed-point iteration scheme. The coupling begins by initializing $T_l$ to air temperature ($T_a$) and setting $C_i$ equal to ambient CO₂ concentration ($C_a$). It then computes $A_n$ and $G_s$ using the photosynthesis-conductance coupling. These outputs are passed to the Monteith energy balance model, which updates $T_l$ along with sensible, latent (i.e. evapotranspiration), and radiative fluxes. The process iterates until the temperature difference between successive steps falls below a user-defined threshold ($\Delta T_l < \varepsilon$).

While the fixed-point method is robust in most cases, it does not guarantee convergence or uniqueness of the solution under all conditions. Users are advised to monitor convergence diagnostics and reject non-physical results when necessary (e.g. when $T_l$ diverges). In practice, convergence is typically fast but can depend on both environmental conditions and parameter values.

This hard-coupling approach ensures that CO₂ and water fluxes are consistent with leaf thermodynamics. It is flexible enough to accommodate different model combinations, for example, replacing Medlyn's stomatal model by the one of Tuzet requires no changes to the overall iteration logic.

## 2.5 Evaluation strategy

**2.5.1 Datasets.** A first dataset from Schymanski and Or (2017) was used to test the correctness of the coupled energy-balance and heat-transfer formulations. It provides simultaneous measurements of net radiation, latent and sensible heat fluxes for a leaf-sized plate (i.e. artificial leaf) in a controlled wind tunnel, with wind speed varying between 0.9 and 5 m s⁻¹. This artificial leaf has no assimilation rate, making this dataset ideal for isolating physical energy-balance processes without confounding physiological responses. By forcing the same conditions and parameters in PlantBiophysics.jl, we could directly compare simulated and measured fluxes, ensuring that the implementation reproduces known experimental measurements before applying it to real leaves with active gas exchange.

A second dataset was used to evaluate the accuracy and computational performance of the fully coupled energy-balance, photosynthesis and stomatal conductance formulations in

PlantBiophysics.jl. This dataset consists in publicly available gas-exchange measurements of *Eucalyptus delegatensis* collected by Medlyn et al. (2015) at the Tumbarumba eddy covariance site in New South Wales, Australia (Leuning et al. 2005). The data were acquired using a LI-COR 6400 portable gas exchange system during three field campaigns (November 2001, February 2002, and May 2002). The dataset includes spot measurements and A-Ci and A-PAR response curves. The measurement conditions span a wide range of photosynthetically active radiation (PAR, 0–1500 µmol m⁻² s⁻¹) and CO₂ concentrations ($C_a$, 0–1200 ppm), allowing a comprehensive test of the coupled models across light and carbon gradients. Note that the dataset does not include explicit $G_s$-VPD response curves, hence $G_s$ was forced for the simulation of spot measurements for validation purposes.

The spot measurements were used to evaluate simulation behaviour under realistic variability. We used PlantBiophysics.jl to simulate six spot measurements of a single tree for 1 day (14 March 2001), incorporating uncertainty propagation on key input variables. The simulated uncertainty range was derived using MonteCarloMeasurements.jl (Carlson 2020) with input uncertainties of 10% normally distributed error on incident radiation, and a uniform sampling of characteristic leaf dimensions between 0.01 and 0.10 m, ±10 m s⁻¹ air flux, and other uncertainties communicated by LI-COR: ±0.1°C air temperature, ±0.1% pressure, ±1% relative humidity.

A third dataset is used to showcase upscaling at the plant level. This dataset consists in a 3D mock-up of an isolated coffee tree (Pérez-Molina 2018) located at 15° latitude, 100 m altitude, on 12 June 2016, at 12:00 a.m.. Air temperature, wind speed, humidity, and CO₂ concentration are set to 25°C, 1 m s⁻¹, 60%, and 380 ppm, respectively.

**2.5.2 Comparison with other tools from the community and benchmarking.** PlantBiophysics.jl is compared with two commonly used packages for leaf-level ecophysiological modelling: plantecophys (Duursma 2015), a popular R package with a widespread use and a strong reputation in the community, and LeafGasExchange.jl (Yun et al. 2020), a Julia package that implements ecophysiological models using Cropbox.jl (Yun and Kim 2022), a domain-specific language for crop models.

All three tools were used to simulate the dataset from Medlyn et al. (2015) with the same coupled system of photosynthesis, stomatal conductance, and energy balance, comparing four key variables: leaf temperature ($T_b$, °C), net carbon assimilation ($A_n$, µmol CO₂ m⁻² s⁻¹), transpiration ($E$, mmol H₂O m⁻² s⁻¹), and stomatal conductance ($G_s$, µmol CO₂ m⁻² s⁻¹). The system was simulated using 'run!' with a 'ModelList' defined as described in Fig. 6 for PlantBiophysics.jl, 'photosynEB' for plantecophys and 'simulate' with 'ModelC3MD' for LeafGasExchange.jl. In all three cases, Medlyn's model was used for stomatal conductance. For better comparability, the unit propagation system in LeafGasExchange.jl was disabled.

To ensure a fair comparison, the photosynthesis parameters (VcMax, JMax, TPU, and Rd) were fitted using PlantBiophysics.jl and the same values are used in all three packages. Common values were also used for all the other parameters. While parameters estimate differed slightly when fitting with 'fitaci'
<table><caption>**Table 2.** Ranges of values considered as representative and used for benchmarking the three packages.</caption><thead><tr><th>Variable</th><th>Description</th><th>Unit</th><th>Minimum</th><th>Maximum</th></tr></thead><tbody><tr><td>R<sub>s</sub></td><td>Short-wave radiation</td><td>W m<sup>-2</sup></td><td>10.0</td><td>500.0</td></tr><tr><td>T</td><td>Air temperature</td><td>°C</td><td>18.0</td><td>40.0</td></tr><tr><td>Wind</td><td>Wind velocity</td><td>m s<sup>-1</sup></td><td>0.5</td><td>20.0</td></tr><tr><td>P</td><td>Air pressure</td><td>kPa</td><td>90.0</td><td>101.0</td></tr><tr><td>Rh</td><td>Relative humidity</td><td>0–1</td><td>0.1</td><td>0.98</td></tr><tr><td>C<sub>a</sub></td><td>CO<sub>2</sub> air concentration</td><td>ppm</td><td>360.0</td><td>900.0</td></tr><tr><td>SkyFraction</td><td>View factor between the object and the sky for both faces</td><td>0–2</td><td>0.0</td><td>2.0</td></tr><tr><td>d</td><td>Characteristic length</td><td>m</td><td>0.001</td><td>0.5</td></tr><tr><td>JMaxRef</td><td>Potential rate of electron transport</td><td>µmol m<sup>-2</sup> s<sup>-1</sup></td><td>200</td><td>300</td></tr><tr><td>VcMaxRef</td><td>Maximum rate of electron transport</td><td>µmol m<sup>-2</sup> s<sup>-1</sup></td><td>150</td><td>250</td></tr><tr><td>RdRef</td><td>Mitochondrial respiration in the light</td><td>µmol m<sup>-2</sup> s<sup>-1</sup></td><td>0.3</td><td>2.0</td></tr><tr><td>TPURef</td><td>Triose-phosphate utilization-limited photosynthesis rate</td><td>µmol m<sup>-2</sup> s<sup>-1</sup></td><td>5.0</td><td>20.0</td></tr><tr><td>g0</td><td>Stomatal conductance parameter to CO<sub>2</sub></td><td>µmol m<sup>-2</sup> s<sup>-1</sup></td><td>0.001</td><td>2.0</td></tr><tr><td>g1</td><td>Stomatal conductance parameter to CO<sub>2</sub></td><td>–</td><td>0.5</td><td>15.0</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Energy balance measured on an artificial leaf (points) and simulated by PlantBiophysics.jl (lines) using the model of Monteith and Unsworth (2013) corrected by Schymanski and Or (2017). Data from Schymanski et al. (2017, Fig. 6A). LE, latent heat flux; H, sensible heat flux; Rn, net radiation flux; H + LE, sum of H and LE components. The assimilation is forced to 0.0 in the simulation using the 'ConstantA' model.</figcaption></figure>

and 'fitBB' in plantecophys, the resulting comparisons were found to be of similar magnitude (not shown).

For computational benchmarking, we generated 100 random input conditions within representative ranges (see Table 2). These included values for air temperature, wind speed, radiation, VPD, CO₂ concentration, and photosynthetic parameters. We then ran 100 simulations per condition using the BenchmarkTools.jl package (Chen and Revels 2016) for Julia tools, and microbenchmark for R. All tests were performed on a MacBook Pro with an Apple M3 Max processor (14 cores) and 36 Go RAM.

Each simulated variable was compared with observations using common performance metrics: the Nash–Sutcliffe model

efficiency coefficient (EF), the root-mean-squared error (RMSE), and its normalized version (nRMSE), and the bias (Bias) and its normalized version (nBias).

## 3. USE CASES

### 3.1 Evaluation with an artificial leaf

The simulation of the artificial leaf reproduced the general magnitude and trends of the measured fluxes (Fig. 7): net radiation remained stable across wind speeds, latent heat flux increased, and sensible heat flux decreased, reflecting the observed enhanced evaporative cooling. Minor discrepancies with
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Measured (x-axis) and simulated (y-axis) *E. delegatensis* (Medlyn et al. 2015) net carbon assimilation (A, µmol m⁻² s⁻¹), transpiration rate (B, mmol m⁻² s⁻¹), stomatal conductance for CO₂ (C, µmol m⁻² s⁻¹) and leaf temperature (D, °C). All simulations were performed using a photosynthesis-stomatal conductance-energy balance coupled model with LeafGasExchange.jl (green circles), plantecophys (blue) and PlantBiophysics.jl (red). Grey line represents *x* = *y*.</figcaption></figure>

observations are consistent with experimental uncertainties and simplifications in leaf property parameterization. Overall, this evaluation indicates that the model implementation is physically consistent and reproduces the underlying processes observed in controlled experiments.

## 3.2 Packages comparison

PlantBiophysics.jl simulated net CO₂ assimilation (A_n) with high accuracy, with an nRMSE of 5% (RMSE: 2.19 µmol m⁻² s⁻¹, EF: 0.96). In comparison, plantecophys produced higher error (nRMSE: 9%, RMSE: 4.03 µmol m⁻² s⁻¹, EF: 0.86) with a skewed prediction trend, while LeafGasExchange.jl simulations had more than twice the error (12%, 5.59 µmol m⁻² s⁻¹, 0.72, respectively) with a tendency

to underestimate values in the 10–30 µmol m⁻² s⁻¹ range (Fig. 8A).

For transpiration rates (E), PlantBiophysics.jl showed similar errors overall (nRMSE: 6%, RMSE: 0.5 mmol m⁻² s⁻¹, EF: 0.89). plantecophys exhibited higher dispersion (8%, 0.73 mmol m⁻² s⁻¹, 0.77, respectively), while LeafGasExchange.jl had a wider spread and lower accuracy (11%, 0.94 mmol m⁻² s⁻¹, 0.62, respectively, Fig. 8B).

Simulations of stomatal conductance to CO₂ (G_s) were satisfactory for PlantBiophysics.jl (RMSE: 0.02 µmol m⁻² s⁻¹, EF: 0.86) while plantecophys showed a larger error (0.08, -0.53, respectively) with a tendency to overestimation (+14% bias), and LeafGasExchange.jl had a very high error (0.05, 0.39, respectively) with an overall underestimation bias of -10% (Figs 8C and 9).
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Taylor diagram (Taylor 2001) illustrating the correlation (patterns similarity) and standard deviation between simulations from PlantBiophysics.jl (red), plantecophys (blue) and LeafGasExchange.jl (green) and data from *E. delegatensis* (Medlyn et al. 2015) for net carbon assimilation (circles, µmol m⁻² s⁻¹), net transpiration rate (squares, mmol m⁻² s⁻¹), leaf temperature (diamond, °C) and stomatal conductance for CO₂ (triangle, µmol m⁻² s⁻¹). Standard deviations are normalized with respect to data (black circle). The closer a point is from the observed data (black circle) the better.</figcaption></figure>

Leaf temperature ($T_l$) was accurately predicted by all three packages. PlantBiophysics.jl achieved the lowest RMSE (0.51°C) and highest EF (0.99), followed by LeafGasExchange.jl (RMSE: 0.72°C, EF: 0.97), while planteco-phys showed slightly higher errors (RMSE: 1.56°C, EF: 0.88, Fig. 8D).

These results are synthesized in a Taylor diagram (Fig. 9), which represents both the correlation and standard deviation of each model's predictions relative to the observations. Across all four variables: net assimilation, transpiration, stomatal conductance, and leaf temperature, PlantBiophysics.jl simulations aligned closest with the observed values, both in terms of pattern similarity and amplitude.

These results confirm that PlantBiophysics.jl reliably simulates coupled biophysical processes of energy balance, photosynthesis, and stomatal conductance. In direct comparison with the two other tools, PlantBiophysics.jl consistently produced lower errors and higher modelling efficiency. This improved performance likely reflects both the incorporation of updated process representations such as the correction to the Penman-Monteith model by Schymanski and Or (2017), and the structured design for model implementation in PlantSimEngine.jl, which reduces the risk of errors during model coupling. A more detailed investigation is required to precisely identify the sources of errors, but this is beyond the scope of this study that is primarily focused on demonstrating that the

implementations of the models provided by the package are extensively tested and produce results consistent with those from other widely used tools in the community.

### 3.3 Uncertainty propagation in simulations

The simulation of the six spot measurements showed good agreement with the observations (Table 3) of leaf-to-air vapour pressure deficit ($D_l$, nRMSE: 8.5%), leaf temperature ($T_l$, nRMSE: 4.6%), assimilation (A, nRMSE: 21.6%) and transpiration (Tr, nRMSE: 4.2%). Assimilation was overestimated at two time points (15:03 and ~16:18), which contributed to an overall positive bias of 15.8% (Fig. 10). The propagated uncertainty bands remained narrow for most variables, with standard deviations below 2% of the mean across simulations, indicating relatively stable model behaviour under input variability.

The capacity for propagating uncertainty throughout the simulation chain offers a powerful feature that is rarely integrated in other ecophysiological modelling tools. This approach enables users to evaluate how measurement variability or uncertain parameters affect key outputs at the scale of the leaf, or at larger scales if integrated in a multiscale design. Altogether, they can support more informed model interpretation or experimental design. While Monte Carlo methods used for uncertainty propagation are computationally intensive, their use is well supported by the performance of the package.
<table><caption>**Table 3.** Model performance statistics related to simulations presented in Fig. 10.</caption><thead><tr><th>Variable</th><th>RMSE</th><th>nRMSE (%)</th><th>EF</th><th>Bias</th><th>nBias (%)</th></tr></thead><tbody><tr><td>D<sub>l</sub> (kPa)</td><td>0.078 ± 0.016</td><td>8.5 ± 1.7</td><td>0.95 ± 0.019</td><td>0.056 ± 0.018</td><td>6.1 ± 2.0</td></tr><tr><td>T<sub>l</sub> (°C)</td><td>0.45 ± 0.074</td><td>4.6 ± 0.85</td><td>0.98 ± 0.0079</td><td>0.33 ± 0.11</td><td>3.5 ± 1.1</td></tr><tr><td>A (µmol m<sup>-2</sup> s<sup>-1</sup>)</td><td>3.15 ± 0.25</td><td>21.6 ± 1.8</td><td>0.68 ± 0.051</td><td>21 ± 0.24</td><td>15.8 ± 1.7</td></tr><tr><td>Tr (mol m<sup>-2</sup> s<sup>-1</sup>)</td><td>0.14 ± 0.046</td><td>4.2 ± 1.4</td><td>0.98 ± 0.012</td><td>083 ± 0.032</td><td>2.5 ± 1.0</td></tr></tbody></table>

RMSE, root mean squared error; nRMSE, normalized RMSE; EF, modelling efficiency; nBias, normalized bias; D<sub>b</sub>, leaf-to-air vapour pressure deficit; T<sub>b</sub>, leaf temperature; A, assimilation; Tr, transpiration.

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Leaf-to-air vapour pressure deficit (D<sub>l</sub>), leaf temperature (T<sub>l</sub>), carbon assimilation (A) and transpiration (Tr) of a leaf in snap measurements from Medlyn et al. (2015). Measurements are denoted in red circles and PlantBiophysics.jl simulations in blue lines ± propagated error from 10% normally distributed error on radiation, 0.1°C air temperature, 10 m s<sup>-1</sup> air flux, 0.1% pressure, 1% relative humidity and 0.01–0.10 uniform sampling of leaf characteristic dimension. The photosynthesis is simulated using the model from Farquhar et al. (1980) with the analytical resolution proposed by Baldocchi (1994) and fitted on A – C<sub>i</sub> curves data. The energy balance is simulated using the model of Monteith and Unsworth (2013) corrected by Schymanski and Or (2017). The stomatal conductance is forced to the measurement value.</figcaption></figure>

Linear uncertainty propagation methods are also available through the Measurements.jl package that greatly reduce computation time, but is generally not suitable for the highly nonlinear nature of biophysical models. Therefore, the choice of uncertainty propagation method should be carefully considered based on the specific characteristics of the models and the desired level of accuracy.

Overall, uncertainty propagation is a particularly interesting and freely available feature made possible by the ability to pass any type as input, including uncertainty-aware numerical types. Ultimately, its use is entirely optional and left up to the user, depending on the context and objectives of the analysis.

## 3.4 Computational performance

Simulations of the full energy balance–photosynthesis–stomatal conductance coupling such as results displayed in Fig. 11 had a

computational mean time of 5.8 µs with standard deviation of 3.8 µs and a median computational time of 5.3 µs for PlantBiophysics.jl. The highest computational time was reached in 77.9 µs while the lower in ~4.9 µs. A simulation on one time step and one leaf did 181 allocations for an estimated 14.41 KiB.

As for the other two packages, plantecophys had a computational mean time of 39.4 ms with standard deviation of 6.0 ms and a median computational time of 38.8 ms. The highest computational time was reached in 154.8 ms and the lower in 30.6 ms. LeafGasExchange.jl had a computational mean time of 2.7 ms with standard deviation of 1.6 ms and a median computational time of 1.9 ms. The highest computational time was reached in 35.7 ms and the lower in 1.7 ms.

On average, a coupled energy balance–photosynthesis–stomatal conductance simulation with PlantBiophysics.jl ran 6739 times faster than the same set-up with plantecophys and 462 times faster than LeafGasExchange.jl (averaging factor between each run).
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Histograms representing the computational time used for the same set of random input conditions with (A) PlantBiophysics.jl, (B) R package plantecophys and (C) LeafGasExchange.jl. Note that x-axis is logarithmic, and y-axis represents the density (i.e. the fraction of simulations). Mean computational times (red dashed lines) and standard deviation intervals (orange background) are also plotted.</figcaption></figure>

This notable performance difference can be attributed to different factors. First, the inherent performance characteristic of the programming languages used plays a significant role. PlantBiophysics.jl is developed in Julia, which provides fast execution speed. In contrast, the plantecophys package is implemented in R, an interpreted language valued for its ease of use and wide range of available packages, but not specifically recognized for its speed. LeafGasExchange.jl, although implemented in Julia, uses Cropbox.jl, which prioritizes providing a simple and standardized approach for model implementation and coupling, rather than emphasizing performance. As a result, it may not enforce strict type stability, which can impact computational efficiency. This shows that albeit Julia being highly performant, users should care about how their models are implemented to ensure type stability, and hence performance.

Simulation time in PlantBiophysics.jl can appear to scale more efficiently with the number of simulated objects than with the number of time steps. For example, simulating 1000 objects for a single time step is roughly four times faster than simulating a single object for 1000 time-steps (Table 4). This difference arises from the fixed initialization cost incurred when multiple time steps are used, due to input checks and memory pre-allocation. However, this overhead is quickly amortized at larger scales: simulating 1000 time-steps over 1000 leaves still takes the same time per unit of simulation as simulating one object over 1000 time-steps.

Additionally, the impressive performance of PlantBiophysics.jl can be attributed to the design principles of PlantSimEngine.jl, which prioritizes performance considerations from the start.

<table><caption>**Table 4.** Simulation duration over time-steps and objects.</caption><thead><tr><th>Object type</th><th>Num. of objects</th><th>Num. of time-steps</th><th>Total time</th><th>Time per unit (1 object, 1 time-step) (μs)</th></tr></thead><tbody><tr><td>Leaf</td><td>1</td><td>1</td><td>5.7 μs</td><td>5.7</td></tr><tr><td>Leaf</td><td>1000</td><td>1</td><td>7.5 ms</td><td>7.5</td></tr><tr><td>Leaf</td><td>1</td><td>1000</td><td>29.2 ms</td><td>29.2</td></tr><tr><td>Leaf</td><td>1000</td><td>1000</td><td>30 s</td><td>30.0</td></tr><tr><td>Plant (leaves)</td><td>2179</td><td>3</td><td>25.5 ms</td><td>3.9</td></tr></tbody></table>

PlantSimEngine.jl offers higher-level methods for managing model coupling using a dependency graph, as well as checks and parallelization over time-steps and objects. These design choices help optimize computational efficiency and contribute to the superior performance observed in PlantBiophysics.jl, with no additional work required from the user.

## 3.5 Example use case: model application on a 3D scene

PlantBiophysics.jl is highly composable and can be used in any type of workflow, including parameter fitting, simulation and visualization. Figure 12 demonstrates a simulation on a 3D mock-up of a coffee tree leveraging the combined capabilities of PlantBiophysics.jl, PlantSimEngine.jl, PlantMeteo.jl, MultiScaleTreeGraph.jl, and PlantGeom.jl. We used PlantBiophysics.jl to simulate the biophysical processes at the leaf level using the coupled energy balance—photosynthesis—
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Leaf temperature (°C) of a 3D coffee plant simulated using PlantBiophysics.jl. The photosynthesis is simulated using the model from Farquhar et al. (1980) with the analytical resolution proposed by Baldocchi (1994). The energy balance is simulated using the model of Monteith and Unsworth (2013) corrected by Schymanski and Or (2017). The stomatal conductance is simulated using the model of Medlyn et al. (2015). The *x*, *y* and *z* axes are given in metre.</figcaption></figure>

stomatal conductance models described earlier. The simulation took 25.5 ± 2.8 ms (47.9 ms on Github's workers), in sequential mode over 2179 leaves and three time-steps, with an average of 3.9 µs per leaf. This improvement in the simulation time per node and per leaf can be explained by the setup of uniform model values between all leaves compared with the more complex setup above where each leaf has its own parameters, gaining time in the initialization step, and because the computations on a plant graph takes a different, more optimized call path.

## 4. CONCLUSION

Many existing tools for simulating biophysical processes strike a trade-off between usability and performance: interactive and user-friendly packages are well suited for model calibration and exploratory simulations but are often too slow for large-scale or long-term simulations; conversely, high-performance tools are typically more complex to use and less flexible for rapid prototyping. This limits the ability to integrate detailed biophysical models into realistic simulation setups involving millions of leaves and numerous time steps.

PlantBiophysics.jl addresses this gap by combining ease of use with high computational performance, leveraging the Julia programming language and the PlantSimEngine.jl modelling framework to deliver several orders of magnitude faster computations

than comparable tools such as plantecophys (R) and LeafGasExchange.jl (Julia). This speed makes it possible to scale from leaf-level measurements to full-plant or plot-scale simulations without compromising interactivity or flexibility. Its composable, type-generic design also allows seamless integration with other packages, enabling advanced features such as uncertainty propagation.

Built on PlantSimEngine.jl, the software architecture supports easy model definition, coupling, and swapping. Models can be swapped or combined with minimal code changes, facilitating intercomparison and targeted testing.

Evaluation with field data confirmed strong performance. It was consistently closer to the observations compared with other packages, with nRMSE below or equal to 7% for leaf-to-air vapour pressure deficit, leaf temperature, assimilation and transpiration. It also provided modelling efficiency between 0.68 and 0.98 when applied on independent time-steps during a day.

By bridging the gap between user-friendliness and computational scalability, PlantBiophysics.jl provides a unified solution for simulating, comparing, and scaling biophysical models. This makes it particularly valuable for research that requires both detailed leaf-level process representation and operational-scale simulation capacity.

Future developments may focus on expanding the model library, better integration with the other packages from the Virtual Plant Laboratory (Morales et al. 2025), coupling it with other tools from the community, and ensuring
compatibility with other advanced modelling frameworks such as OpenAlea (Pradal et al. 2008) or ModelingToolkit.jl (Ma et al. 2021). The package will also be used in broader models, for example to account for the effect of microclimate on plant growth and development in complex systems.

## AUTHOR CONTRIBUTIONS

Rémi Vezy (Conceptualization, Methodology, Software, Validation, Formal analysis, Investigation, Data curation, Writing—original draft, Writing—review & editing, Visualization, Supervision, Project administration, Resources, Funding acquisition), Simon Treillou (Methodology, Software, Validation, Formal analysis, Investigation, Visualization, Writing—original draft, Writing—review & editing), Samuel Mackeown (Software, Resources), Emilie Peynaud (Supervision, Writing—review & editing), Jean Dauzat (Resources, Funding acquisition), Raphaël P.A. Perez (Resources, Funding acquisition), and Thomas Arsouze (Writing—review & editing, Resources)

Conflict of interest: None declared.

## FUNDING

This work was supported by the AMAP laboratory 'Botany and Modelling of Plant Architecture and Vegetation', which funded Simon Treillou's internship in 2021. It also received partial support from the PalmStudio research project, co-funded by the SMART research institute and CIRAD (French Agricultural Research Centre for International Development). The authors have been supported by the MaCS4Plants CIRAD network, initiated from the AGAP Institute and AMAP joint research units.

## DATA AVAILABILITY

This study used PlantBiophysics.jl 0.16.2, available from the General Julia registry, its Github repository (https://github.com/VEZY/PlantBiophysics.jl), and archived on Zenodo (Vezy et al. 2025). Its documentation is available online at https://vezy.github.io/PlantBiophysics.jl. All results are fully reproducible using the scripts and notebooks provided in https://github.com/VEZY/PlantBiophysics-paper, which are automatically rendered at https://vezy.github.io/PlantBiophysics-paper. Developers of other tools are welcome to contribute comparisons by submitting a pull request to the repository.

## REFERENCES

* Bailey BN, Ponce de León MA, Krayenhoff ES. One-dimensional models of radiation transfer in heterogeneous canopies: a review, re-evaluation, and improved model. *Geosci Model Dev* 2020;13:4789–808. https://doi.org/10.5194/gmd-13-4789-2020

* Baldocchi D. An analytical solution for coupled leaf photosynthesis and stomatal conductance models. *Tree Physiol* 1994;14:1069–79. https://doi.org/10.1093/treephys/14.7.8-9.1069

* Bezanson J, Edelman A, Karpinski S, et al. Julia: a fresh approach to numerical computing. *SIAM Rev* 2017;59:65–98. https://doi.org/10.1137/141000671

* Brisson N, Mary B, Ripoche D, et al. STICS: a generic model for the simulation of crops and their water and nitrogen balances. I. Theory and parameterization applied to wheat and corn. *Agronomie* 1998;18:311–46. https://doi.org/10.1051/agro:19980501

* Buckley TN. Modeling stomatal conductance. *Plant Physiol* 2017;174:572–82. https://doi.org/10.1104/pp.16.01772

* Carlson FB. MonteCarloMeasurements.jl: nonlinear propagation of arbitrary multivariate distributions by means of method overloading. *arXiv*, https://doi.org/10.48550/arXiv.2001.07625, 21 January 2020, preprint: not peer reviewed.

* Chen J, Revels J. Robust benchmarking in noisy environments. *arXiv*, https://doi.org/10.48550/arXiv.1608.04295, 15 August 2016, preprint: not peer reviewed.

* Crous KY. Plant responses to climate warming: physiological adjustments and implications for plant functioning in a future, warmer world. *Am J Bot* 2019;106:1049–51. https://doi.org/10.1002/ajb.1329

* Damour G, Simonneau T, Cochard H, et al. An overview of models of stomatal conductance at the leaf level. *Plant Cell Environ* 2010;33:1419–38. https://doi.org/10.1111/j.1365-3040.2010.02181.x

* Duursma RA. Plantecophys—an R package for analysing and modelling leaf gas exchange data. *PLoS One* 2015;10:e0143346. https://doi.org/10.1371/journal.pone.0143346

* Duursma RA, Medlyn BE. MAESPA: a model to study interactions between water limitation, environmental drivers and vegetation function at tree and stand levels, with an example application to [CO₂] × drought interactions. *Geosci Model Dev* 2012;5:919–40. https://doi.org/10.5194/gmd-5-919-2012

* Ely KS, Rogers A, Agarwal DA, et al. A reporting format for leaf-level gas exchange data and metadata. *Ecol Inform* 2021;61:101232. https://doi.org/10.1016/j.ecoinf.2021.101232

* Evers JB, Letort V, Renton M, et al. Computational botany: advancing plant science through functional-structural plant modelling. *Ann Bot* 2018;121:767–72. https://doi.org/10.1093/aob/mcy050%

* Farquhar GD, von Caemmerer S, Berry JA. A biochemical model of photosynthetic CO₂ assimilation in leaves of C₃ species. *Planta* 1980;149:78–90. https://doi.org/10.1007/BF00386231

* Giordano M. Uncertainty propagation with functionally correlated quantities. *arXiv*, https://doi.org/10.48550/arXiv.1610.08716, 27 October 2016, preprint: not peer reviewed.

* Godin C, Caraglio Y. A multiscale model of plant topological structures. *J Theor Biol* 1998;191:1–46. https://doi.org/10.1006/jtbi.1997.0561

* Gutschick VP. Leaf energy balance: basics, and modeling from leaves to canopies. In: Hikosaka K, Niinemets Ü, Anten NPR (eds.), *Canopy Photosynthesis: From Basics to Applications, Advances in Photosynthesis and Respiration*. Dordrecht, Netherlands: Springer, 2016, 23–58. https://doi.org/10.1007/978-94-017-7291-4_2

* Leuning R, Cleugh HA, Zegelin SJ, et al. Carbon and water fluxes over a temperate Eucalyptus forest and a tropical wet/dry savanna in Australia: measurements and comparison with MODIS remote sensing estimates. *Agric For Meteorol* 2005;129:151–73. https://doi.org/10.1016/j.agrformet.2004.12.004

* Liu Y, Wang X, Fan D, et al. The use of R in photosynthesis research. *Funct Plant Biol* 2022;49:565–72. https://doi.org/10.1071/FP21102

* Lombardozzi DL, Bonan GB, Smith NG, et al. Temperature acclimation of photosynthesis and respiration: a key uncertainty in the carbon cycle-climate feedback. *Geophys Res Lett* 2015;42:8624–31. https://doi.org/10.1002/2015GL065934

* Lombardozzi DL, Smith NG, Cheng SJ, et al. Triose phosphate limitation in photosynthesis models reduces leaf photosynthesis and global terrestrial carbon storage. *Environ Res Lett* 2018;13:074025. https://doi.org/10.1088/1748-9326/aacf68

* Ma Y, Gowda S, Anantharaman R, et al. ModelingToolkit: a composable graph transformation system for equation-based modeling. *arXiv*, https://doi.org/10.48550/arXiv.2103.05244, 9 March 2021, preprint: not peer reviewed.
Medlyn B, Pepper D, Keith H. Tumbarumba gas exchange. 2015. https://doi.org/10.6084/m9 figshare.1538079.v1

Medlyn BE, Duursma RA, Eamus D, et al. Reconciling the optimal and empirical approaches to modelling stomatal conductance. *Glob Chang Biol* 2011;17:2134–44. https://doi.org/10.1111/j.1365-2486.2010.02375.x

Monteith JL, Unsworth MH. Chapter 13-steady-state heat balance:(i) water surfaces, soil, and vegetation. In: Monteith JL, Unsworth MH (eds.), *Principles of Environmental Physics*, 4th edn. Boston: Academic Press, 2013, 217–47.

Morales A, Kottelenberg DB, Ernst A, et al. The Virtual Plant Laboratory: a modern plant modeling framework in Julia. *In Silico Plants* 2025;7:diaf005. https://doi.org/10.1093/insilicoplants/diaf005

Muir CD. Tealeaves: an R package for modelling leaf temperature using energy budgets. *AoB Plants* 2019;11:plz054. https://doi.org/10.1093/aobpla/plz054

Parkhurst DF, Loucks OL. Optimal leaf size in relation to environment. *J Ecol* 1972;60:505–37. https://doi.org/10.2307/2258359

Pérez-Molina JP. Impacts of drought on coffee: integrating physiological and morphological processes from the leaf to the whole-plant scale. Universidade Federal de Viçosa. 2018.

Ponce de León MA, Bailey BN. Evaluating the use of Beer's law for estimating light interception in canopy architectures with varying heterogeneity and anisotropy. *Ecol Modell* 2019;406:133–43. https://doi.org/10.1016/j.ecolmodel.2019.04.010

Pradal C, Dufour-Kowalski S, Boudon F, et al. OpenAlea: a visual programming and component-based software platform for plant modeling. *Funct Plant Biol* 2008;35:751–60. https://doi.org/10.1071/FP08084

Riley WJ, Zhu Q, Tang JY. Weaker land-climate feedbacks from nutrient uptake during photosynthesis-inactive periods. *Nat Clim Change* 2018;8:1002–6. https://doi.org/10.1038/s41558-018-0325-4

Roesch E, Greener JG, MacLean AL, et al. Julia for biologists. *Nat Methods* 2023;20:655–64. https://doi.org/10.1038/s41592-023-01832-z

Rogers A, Medlyn BE, Dukes JS, et al. A roadmap for improving the representation of photosynthesis in Earth system models. *New Phytol* 2017;213:22–42. https://doi.org/10.1111/nph.14283

Schymanski SJ, Or D. Leaf-scale experiments reveal an important omission in the Penman-Monteith equation. *Hydrol Earth Syst Sci* 2017;21:685–706. https://doi.org/10.5194/hess-21-685-2017

Stinziano JR, Roback C, Sargent D, et al. Principles of resilient coding for plant ecophysiologists. *AoB Plants* 2021;13:plab059. https://doi.org/10.1093/aobpla/plab059

Taylor KE. Summarizing multiple aspects of model performance in a single diagram. *J Geophys Res* 2001;106:7183–92. https://doi.org/10.1029/2000JD900719

Tuzet A, Perrier A, Leuning R. A coupled model of stomatal conductance, photosynthesis and transpiration. *Plant Cell Environ* 2003;26:1097–116. https://doi.org/10.1046/j.1365-3040.2003.01035.x

Vezy R, le Maire G, Christina M, et al. DynACof: a process-based model to study growth, yield and ecosystem services of coffee agroforestry systems. *Environ Model Softw* 2020;124:104609. https://doi.org/10.1016/j.envsoft.2019.104609

Vezy R. MultiScaleTreeGraph.jl: read, write, analyze, compute and plot multi-scale tree graph files. Version v0.14.4. [Computer software]. Zenodo. https://doi.org/10.5281/zenodo.15832373. 7 July 2025.

Vezy R. PlantSimEngine: a simulation engine for the soil-plant-atmosphere system. *J Open Source Softw* 2023;8:5371. https://doi.org/10.21105/joss.05371

Vezy R, Treillou S, Mackeown S, et al. VEZY/PlantBiophysics.jl. Version v0.16.2. [Computer software]. Zenodo. https://doi.org/10.5281/zenodo.16848384. 13 August 2025.

von Caemmerer S, Farquhar GD. Some relationships between the biochemistry of photosynthesis and the gas exchange of leaves. *Planta* 1981;153:376–87. https://doi.org/10.1007/BF00384257

Wang Y, Köhler P, He L, et al. Testing stomatal models at the stand level in deciduous angiosperm and evergreen gymnosperm forests using CliMA Land (v0.1). *Geosci Model Dev* 2021;14:6741–63. https://doi.org/10.5194/gmd-14-6741-2021

Yun K, Kim S-H. Cropbox: a declarative crop modelling framework. *In Silico Plants* 2022;5:diac021. https://doi.org/10.1093/insilicoplants/diac021

Yun K, Timlin D, Kim S-H. Coupled gas-exchange model for C4 leaves comparing stomatal conductance models. *Plants* 2020;9:1358. https://doi.org/10.3390/plants9101358