Canopy & light-environment models

implementation

Prerequisites

This page builds on The big picture and the FF16 model. Read those first if the notation here is unfamiliar.

The idea first

Imagine looking down on a forest. Tall trees intercept light first, leaving shorter trees in their shade. A leaf lower in the canopy has more leaf area above it, so less open-sky light reaches it. This accumulation of shade with depth is the central idea behind a canopy, or light-environment, model.

Rather than following every leaf, the model constructs a vertical light profile: at every height, it records the fraction of open-sky light that remains. Near the canopy top, that fraction is close to one; lower down, crowns above have filtered more of the light. All plants in a patch use the same profile. Each plant reads it to determine the light reaching its leaves, while also changing it by shading the plants below. Growth therefore depends not only on a plant’s own size and strategy, but also on the light environment created by its neighbours.

plant provides several ways to represent this process. Some resolve shading throughout each crown; others use faster approximations that sacrifice some detail. This page explains the six implemented canopy and shading models, how each calculates the light environment and assimilation, the trade-off between realism and speed, and why the perfect-plasticity approximation (PPA) must be smoothed for plant’s adaptive solver. It is an implementation reference for the current code.

How canopy shading affects the model

Light interception, and why it matters

Light is the currency of forest competition. A plant’s carbon gain depends on the light intercepted by its leaves. That light has already been filtered by taller neighbours’ crowns and by leaves higher in the plant’s own crown (self-shading). plant represents this with one shared patch light profile, $E(z)$: the fraction of open-sky light reaching height $z$. Every plant both reads this profile to photosynthesise and contributes to it, because its leaves lower $E$ for plants beneath them.

This profile is consequential, not a minor implementation detail. It determines plant growth, which determines stand size structure and competitive hierarchy, and ultimately affects which strategies persist and which traits selection favours. Errors in shading can therefore propagate to the ecological and evolutionary predictions that the model is designed to make.

Realism versus speed

A fully detailed model would resolve the light field throughout every crown and integrate photosynthesis across it. This is costly, particularly when leaf-level photosynthesis is itself expensive: TF24 uses coupled stomatal/hydraulic optimisation rather than a simple light-response curve. A common alternative is to simplify the crown to a point, or the canopy to a small number of layers. That saves computation but reduces fidelity.

plant lets you choose a point on this spectrum. The important question is not which simplification is universally best, but what each one changes. As the model comparison shows, an approximation can be suitable for one purpose and seriously misleading for another.

The models

The FF16 strategy offers six shading models, selected with control$shading_model. Start by separating two ideas: the shade a plant casts (competition) and how it turns the light it receives into assimilation. Four models use the same smooth Yokozawa leaf-area profile for competition. They differ only in how they calculate a plant’s own assimilation, except that PPA also changes how the light profile is assembled:

"deep-crown" (the default) – the original plant behaviour. Leaf area is distributed continuously through the crown according to the Yokozawa (1995) foliage profile. Gross assimilation is calculated by integrating photosynthesis over crown depth against the light profile. The resulting patch light profile is smooth with height.
"crown-centre" – builds the light profile exactly as "deep-crown" does, but approximates assimilation with a single evaluation at the crown centre rather than an integral over crown depth. Skipping the per-plant quadrature makes it faster, but removes within-crown self-shading.
"mean-light" – sits between deep-crown and crown-centre. It integrates the light over crown depth, weighted by leaf-area density, to obtain one leaf-area-weighted mean light level, then evaluates photosynthesis once at that mean. It therefore retains the changing within-crown light environment as the canopy develops, while replacing many photosynthesis calculations with one per plant. The saving is modest for FF16, whose leaf photosynthesis is a cheap saturating function. It matters more for models with expensive leaf submodels: TF24’s coupled stomatal/hydraulic optimisation makes integrating light and optimising once the natural choice. The trade-off is bias. Because photosynthesis saturates with light, the rate evaluated at mean light exceeds the depth-averaged rate (Jensen’s inequality). "mean-light" therefore always gives at least as much assimilation as "deep-crown"; they agree only when light is uniform within the crown.
"ppa" – the perfect-plasticity approximation used by large models such as FATES. It divides the patch light profile into discrete canopy layers: light decreases by a fixed optical depth (control$ppa_layer_optical_depth) in each layer. Assimilation is evaluated at the crown centre, as in "crown-centre", but reads this stepped profile. A hard step is not differentiable and causes the ODE solver to fail, so layer boundaries are smoothed with a $C^1$-continuous ramp whose width is set by control$ppa_layer_smoothing. This creates two versions of the single "ppa" model, paralleling the Flat Top pair below: PPA (hard step) (ppa_layer_smoothing = 0) is the literal field discretisation (e.g. Strigul 2008, FATES). It is discontinuous and does not run in plant’s adaptive solver. PPA (smoothed) (ppa_layer_smoothing > 0, default 0.3) is the runnable version.

The remaining two models change competition itself: the shade a plant casts. Instead of spreading leaf area through the crown, they collapse it towards a thin layer at the crown centre. Many large vegetation models, including LPJ-GUESS, use a similar discrete-layer view of canopy competition, so these are useful approximations to examine rather than dismiss:

"flat-top-box" – casts a hard step of shade: full below the crown centre and none above. Across cohorts this makes the patch light profile discontinuous, so plant cannot solve the stand at all.
"flat-top-soft-box" – uses the same box competition but smooths the step into a continuous drop, allowing the light environment to be built. It runs, but its incorrectly shaped shade biases predictions.

The selected model is resolved once when a strategy is prepared, so offering the choice has no per-call cost. TF24 supports "deep-crown", "mean-light" (its default), and "crown-centre", but not the PPA or box variants.

library(plant)
library(dplyr)
library(ggplot2)
library(patchwork)

A patch simulated under each model

We grow one FF16 species in a patch and run the identical patch under each of the four smooth-competition models. Because everything except control$shading_model is fixed, differences among runs can be attributed to the shading model alone.

params <- scm_base_parameters("FF16")
patch  <- add_strategies(params, trait_matrix(0.0825, "lma"), birth_rate = 1)

models <- c("deep-crown", "mean-light", "crown-centre", "ppa")

run_model <- function(model) {
  ctrl <- Control()
  ctrl$shading_model <- model
  run_scm(patch, ctrl = ctrl, collect = TRUE)
}

results <- lapply(models, run_model)
names(results) <- models

Each run returns its full collected state (run_scm(collect = TRUE)): a tidy species table, with one row for every node at every step, and the light environment (env$light_availability). We label each species table with its shading model, then combine the tables so that plant’s tidy helpers can work across all four runs.

species_all <- dplyr::bind_rows(
  lapply(models, function(m) dplyr::mutate(results[[m]]$species, method = m))
)
species_all$method <- factor(species_all$method, levels = models)

# Light-extinction coefficient, used throughout to convert canopy openness to
# cumulative leaf area: E = exp(-k_I * LAI), so LAI_above = -log(E) / k_I.
k_I <- results[[1]]$p$strategies[[1]]$pars$k_I

What happens to the stand

The first and least demanding test is whether a stand develops sensibly. We inspect the same patch in three complementary ways as it grows and self-thins: its size distribution, leaf area index, and vertical light profile.

Size distribution

The size distribution consists of node height trajectories. Each line follows a cohort introduced at a particular patch age; line opacity is scaled by cohort density, so sparse cohorts appear faint. plant’s plot_size_distribution() builds this view directly from the tidy species table, and facets by method make the shading models easy to compare.

p_size <- plot_size_distribution(species_all) +
  facet_wrap(~method) +
  theme(legend.position = "bottom")
p_size

Figure 1: Size-distribution trajectories under each shading model.

The discrete layers in "ppa" create visibly different trajectories. Cohorts grow under nearly constant light within a layer, then accelerate as they emerge into a brighter one.

Leaf area index

Leaf area index (LAI) is total one-sided leaf area per unit ground area. It is the cumulative leaf area above the ground, and the recorded light profile already contains it. By Beer’s law, $E = \exp(-k_I\,\mathrm{LAI})$, so patch LAI is $-\log E(0) / k_I$, where $E(0)$ is canopy openness at ground level. We can therefore read it directly from env$light_availability, without integrating over the size distribution. LAI is real leaf area, so the smooth recorded profile gives its true value for all four models: PPA layering changes the light plants experience, not the leaf area that is present.

lai_all <- dplyr::bind_rows(
  lapply(models, function(m) {
    results[[m]]$env$light_availability |>
      dplyr::group_by(time) |>
      dplyr::slice_min(height, n = 1, with_ties = FALSE) |>
      dplyr::ungroup() |>
      dplyr::transmute(method = m, time, lai = -log(light_availability) / k_I)
  })
)
lai_all$method <- factor(lai_all$method, levels = models)

p_lai <- ggplot(lai_all, aes(time, lai, colour = method)) +
  geom_line(linewidth = 1) +
  labs(x = "Patch age (years)", y = "Leaf area index", colour = "Shading model") +
  theme_classic()
p_lai

Figure 2: Leaf area index over patch age under each shading model.

"deep-crown", "mean-light", and "crown-centre" closely track one another and settle to a steady LAI. Under a nearly uniform light profile, their three ways of summarising within-crown light produce similar carbon gain. "mean-light" and "crown-centre" are slightly richer because they do not fully resolve self-shading. In contrast, "ppa" has large, persistent oscillations: discrete light layers create sharp thresholds for cohort success, leading to successive waves of recruitment and self-thinning. The pulses of emerging cohorts are also visible in the size-distribution panel above.

Vertical light profile

Next, we look within the canopy at selected patch ages. The vertical profile of cumulative leaf area above each height is $\mathrm{LAI}_{above}(z) = -\log E(z) / k_I$, where $E$ is canopy openness. This is the quantity that determines how strongly each individual is shaded.

env$light_availability records the smooth underlying profile shared by all models. PPA plants actually experience a stepped version of that profile, so the "ppa" panel applies the same internal layering transform (see FF16_Environment::step_light).

d <- Control()$ppa_layer_optical_depth   # layer thickness (optical depth)
w <- Control()$ppa_layer_smoothing       # boundary smoothing fraction

# C1-smoothed staircase: flat over the lower (1 - w) of each layer, cubic ramp
# over the top w. Mirrors FF16_Environment::smooth_floor() / step_light().
smooth_floor <- function(u, w) {
  n <- floor(u)
  f <- u - n
  t <- pmax(0, f - (1 - w)) / w
  n + ifelse(f <= 1 - w, 0, t^2 * (3 - 2 * t))
}
ppa_step_light <- function(E, d, w) {
  ifelse(E >= 1 | E <= 0, E,
         exp(-d * smooth_floor(-log(E) / d, w)))
}

target_ages <- c(5, 10, 20, 60)

profile_all <- dplyr::bind_rows(
  lapply(models, function(m) {
    light <- results[[m]]$env$light_availability
    if (m == "ppa") {
      light$light_availability <- ppa_step_light(light$light_availability, d, w)
    }
    # nearest recorded step to each target age
    steps <- results[[m]]$steps
    keep <- sapply(target_ages, function(a) steps$step[which.min(abs(steps$time - a))])
    light |>
      dplyr::filter(step %in% keep) |>
      dplyr::mutate(method = m,
                    age = target_ages[match(step, keep)],
                    lai_above = -log(light_availability) / k_I)
  })
)
profile_all$method <- factor(profile_all$method, levels = models)

p_profile <- ggplot(profile_all,
                    aes(lai_above, height, colour = factor(age), group = age)) +
  geom_path() +
  facet_wrap(~method) +
  labs(x = expression(LAI[above]~(m^2/m^2)), y = "Height (m)",
       colour = "Patch age\n(years)") +
  theme_classic()
p_profile

Figure 3: Vertical profile of cumulative leaf area above each height, at several patch ages.

"deep-crown", "mean-light", and "crown-centre" have the same smooth profile because they construct the light environment identically and differ only in their assimilation calculation. "ppa" has its characteristic staircase: cumulative leaf area is constant within a layer and jumps at layer boundaries.

Computational cost: work per step and solver steps

Speed is why the simpler models exist. One patch is inexpensive, but plant is often run thousands of times for trait sweeps, model fitting, or fitness landscapes. At that scale, reducing the cost of an inner calculation matters.

The saving depends on two competing effects. Unlike many vegetation models, plant integrates a patch with an adaptive ODE solver. The solver chooses large steps where dynamics are smooth and small steps where they are not, aiming to achieve the requested accuracy efficiently. A model’s total cost therefore depends both on the work per step and on the number of steps the solver chooses:

Work per step: "deep-crown" is most expensive because it performs a Gauss-Kronrod quadrature over crown depth for every individual. "crown-centre" and "ppa" replace this with one light evaluation.
Number of steps: a smooth light profile lets the solver stride. PPA’s near-stepped profile is much stiffer, forcing many more and smaller steps.

Which effect dominates depends on how the system is stepped. We therefore report both adaptive cost (the cost normally paid) and per-step cost on a fixed schedule.

# Median wall-clock over a few repeats. `times = NULL` uses the adaptive solver;
# supplying a fixed schedule pins every model to the same steps.
time_model <- function(model, times = NULL, reps = 5) {
  p <- patch
  use_fixed <- !is.null(times)
  if (use_fixed) p$ode_times <- times
  ctrl <- Control()
  ctrl$shading_model <- model
  run1 <- function() run_scm(p, ctrl = ctrl, use_ode_times = use_fixed)
  run1()                                   # warm up
  median(replicate(reps, system.time(run1())[["elapsed"]]))
}

speed_table <- function(timings) {
  data.frame(
    model = names(timings),
    `median time (s)` = round(unname(timings), 4),
    `relative to deep-crown` = round(unname(timings / timings[["deep-crown"]]), 2),
    check.names = FALSE,
    row.names = NULL
  )
}

Adaptive solver: real-world cost

This is plant’s default mode and the cost normally paid. Adaptive stepping is part of what makes plant efficient: the solver spends steps only where the dynamics need them. It also means that cost is not determined solely by work per step. A model that is harder to integrate pays an additional cost through extra steps. Most large vegetation models instead use a fixed timestep, where this second consideration does not arise; we return to that case below.

timings_adaptive <- vapply(models, time_model, numeric(1))
knitr::kable(speed_table(timings_adaptive))

model	median time (s)	relative to deep-crown
deep-crown	0.090	1.00
mean-light	0.089	0.99
crown-centre	0.074	0.82
ppa	0.369	4.10

"crown-centre" is fastest. It shares deep-crown’s smooth light profile, and therefore its solver-step count, but skips the per-plant crown-depth integral. "mean-light" costs about as much as "deep-crown": it still integrates over crown depth, although it integrates light rather than rate, so both per-step work and step count are similar. Despite cheaper per-plant assimilation, "ppa" is slower overall because its near-stepped profile makes the adaptive solver take many more steps. Finer layers (smaller ppa_layer_optical_depth) or sharper boundaries (smaller ppa_layer_smoothing) make this effect stronger.

Fixed ODE schedule: per-step cost

When all models are pinned to the same fixed schedule (run_scm(..., use_ode_times = TRUE)), the step-count effect disappears and only per-step cost remains. PPA’s advantage is then clear: like "crown-centre", it skips crown-depth quadrature, so both are markedly faster than "deep-crown".

times <- seq(0, params$max_patch_lifetime, length.out = 4000)
timings_fixed <- vapply(models, time_model, numeric(1), times = times)
knitr::kable(speed_table(timings_fixed))

model	median time (s)	relative to deep-crown
deep-crown	2.304	1.00
mean-light	2.339	1.02
crown-centre	1.707	0.74
ppa	1.847	0.80

PPA is therefore not inherently slow: on a fixed schedule it is faster than "deep-crown". Its apparent adaptive-solver cost comes entirely from the extra steps required by its stiff profile. This also explains why discrete-layer schemes are unproblematic in large models that use fixed timesteps. There, a stepped, stiff profile costs no additional steps, so the difficulty does not appear. plant’s adaptive solver makes stiffness visible and, as the next section shows, makes a hard step fail outright. "deep-crown" remains the default because it is the most faithful representation of within-crown light capture and has no speed penalty relative to earlier versions of plant.

Why PPA must be smoothed

A literal PPA discretisation is a hard step function. When a growing plant’s crown centre crosses a layer boundary, the light it receives jumps discontinuously, and so does its growth rate. ppa_layer_smoothing exists to remove this problem. The issue is specifically adaptive ODE stepping, not PPA itself.

Why adaptive stepping fails at a discontinuity

plant solves patch dynamics with an adaptive, error-controlled embedded Runge-Kutta scheme. At each step, it makes two state estimates of different orders. Their difference estimates local truncation error. If that error is larger than the requested tolerance, the solver rejects the step and retries with a smaller one. This relies on a smooth solution: halving the step should roughly quarter the error.

That assumption fails at a discontinuity. However small the step becomes, a step that straddles the jump has an O(jump) error that never falls below the tolerance. The solver keeps rejecting and shrinking the step until it gives up. This is exactly what happens with a hard step (ppa_layer_smoothing = 0):

ctrl_hard <- Control()
ctrl_hard$shading_model <- "ppa"
ctrl_hard$ppa_layer_smoothing <- 0      # hard step, no smoothing

tryCatch(
  run_scm(patch, ctrl = ctrl_hard),
  error = function(e) conditionMessage(e)
)

[1] "Non-finite cohort density in the SCM size-density (characteristic) equations: species 1 has a node with density=942350364544097030078472872299071659326557317447943132471173833753559581065216.000000 (log_density=179.542259, height=5.449201) at time=6.788465. The density derivative -d(growth)/d(height) - mortality can grow without bound when growth rate falls steeply with size under rapidly changing or extreme environmental forcing, driving density to overflow. Try a shorter max_patch_lifetime or less extreme environmental drivers."

A fixed schedule (run_scm(..., use_ode_times = TRUE)) avoids the error-control loop because it simply takes its prescribed steps. However, it then has no protection against a large explicit step overshooting the jump, which sends growth to non-finite values unless the schedule is very fine. A hard step is therefore not robustly handled by either approach.

Changing the rounding direction (using light at the top rather than the bottom of a layer) does not solve this. The problem is the discontinuity, not its location.

The smoothed profile

ppa_layer_smoothing replaces each hard boundary with a $C^1$-continuous cubic ramp over the upper fraction of the layer (see FF16_Environment::smooth_floor). The light profile, and therefore each plant’s growth rate, now has a continuous first derivative. The adaptive error estimate behaves as the solver expects, so steps are accepted normally. Even a small amount of smoothing restores a well-behaved integration; with the default, adaptive and fixed-step solutions agree.

Two parameters control the discretisation:

ppa_layer_optical_depth – the optical-depth thickness of a layer (default 0.5, meaning one unit of leaf area index per layer at the default extinction coefficient k_I = 0.5). Smaller values produce finer layers.
ppa_layer_smoothing – the fraction of each layer over which its boundary is smoothed (default 0.3). Values approaching 0 recover a hard step and its instability; 1 removes the flat region and approaches a smooth "deep-crown"-style profile. This creates a trade-off: smaller values more closely represent a true PPA step but are stiffer and slower, as the speed table shows; larger values integrate faster but blur the layers.

How the models represent competition

The four smooth-competition models above construct the same light environment. They differ only in how a plant summarises its own assimilation. The two box models instead alter the shade a plant casts, concentrating leaf area near a thin layer at the crown centre. Because this mirrors the discrete-layer competition used by large vegetation models such as LPJ-GUESS, it is useful to inspect the difference directly.

comp_profile <- function(model, height = 10) {
  s <- FF16_Strategy()
  s$control$shading_model <- model
  ind <- FF16_Individual(s)
  ind$set_state("height", height)
  z <- seq(0, height, length.out = 400)
  data.frame(model = model, z = z,
             competition = sapply(z, function(zz) ind$compute_competition(zz)))
}
comp <- dplyr::bind_rows(
  comp_profile("crown-centre"),
  comp_profile("flat-top-box"),
  comp_profile("flat-top-soft-box")
)

ggplot(comp, aes(competition, z, colour = model)) +
  geom_path() +
  labs(x = "Shade cast (competition) by a 10 m plant", y = "Height (m)",
       colour = "Model") +
  theme_classic()

Figure 4: Shade cast (competition profile) by a 10 m plant under the correct crown-centre profile versus the hard and soft box approximations.

The correct crown-centre profile tapers continuously to zero through the crown. The box profile is a step; the soft box replaces that step with a continuous-but-incorrect drop concentrated near the crown centre.

The hard step does not even run. When summed across a population of discrete cohorts, step competition produces a discontinuous, piecewise-constant patch light profile. The adaptive light-environment spline cannot represent it, so the model fails before a stand develops:

ctrl_box <- Control()
ctrl_box$shading_model <- "flat-top-box"
tryCatch(
  run_scm(patch, ctrl = ctrl_box),
  error = function(e) conditionMessage(e)
)

[1] "Interpolated function as refined as currently possible"

The soft box runs because its continuous profile allows the light environment to be built. However, it distributes shade incorrectly through the canopy and so biases predictions relative to the deep-crown reference.

What this comparison shows

The full model comparison (see the reproduction below) tests these models in increasingly demanding settings: stand dynamics, demographic equilibrium, invasion fitness, and the evolutionary endpoint. Its recurring result is that cheap models are often adequate for stand dynamics and acceptable for equilibrium, whereas invasion fitness is much less forgiving. The key requirement is continuity of the shade each plant casts as a function of height. "deep-crown", "mean-light", and "crown-centre" all construct a continuous light profile and behave well; PPA’s stepped profile and the box models’ incorrectly shaped competition do not.

This does not mean that large models are wrong to use discrete layers. For stand-level simulation on a fixed timestep, that representation is fast, robust, and adequate. The caution applies specifically to the eco-evolutionary questions that plant is built to answer.

The full model comparison and evolutionary analysis – demographic equilibrium, invasion-fitness landscapes, and singular-strategy / work-precision results – is reproduced in the DunFalster 2026 canopy reproduction post.