Self-thinning

example

After disturbance, many young plants may establish together. As they grow, the canopy closes and light becomes unevenly distributed. Taller plants retain access to light, while shaded plants grow slowly and are more likely to die. The stand therefore shifts from many small individuals to fewer large ones.

This decline in density as average size increases is called self-thinning. On logarithmic axes, many stands follow an approximately straight trajectory known as the self-thinning line. Its universality across species, environments, and stand types remains debated.

Here we ask whether this pattern emerges from an FF16 simulation. It is not imposed on the model: we calculate it from the growth and mortality represented by the size-structured PDE.

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

plant_log_console()

We begin with one FF16 strategy and rapidly reduce recruitment after establishment, isolating the development of the initial cohort. Schedule refinement gives the fast early dynamics more resolution:

p0 <- scm_base_parameters("FF16")

p0$strategy_default$pars$recruitment_decay <- 5
# rapidly reduces later recruitment to isolate the first group of plants

ctrl <- Control()
ctrl$schedule_eps <- 0.005
# more refined scheduling (schedule_eps is a Control field, not a Parameters one)

lma_trait <- trait_matrix(0.1457419, "lma")

p1 <- add_strategies(p0, lma_trait)

stand_results <- run_scm(p1, ctrl = ctrl, refine_schedule = TRUE)$parameters
# refines the node-introduction schedule where the estimated error is largest

out <- run_scm(stand_results, ctrl = ctrl, collect = TRUE)
# run the SCM model

The collected species table contains one row per node and time. Expanding its states and integrating over height produces one row of patch totals per time. Dividing size totals by total density then gives mean size per individual:

tidy_stand <- out

stand_expand <- FF16_expand_state(tidy_stand)
# adds in more state variables from our stand output (leaf, bark, sapwood, stem properties to name a few)

patch_total <-
    stand_expand$species %>%
    integrate_over_size_distribution() %>%
    mutate(
      # create average sizes
      across(c(diameter_stem, area_stem, area_leaf, height, mass_above_ground), ~.x/density,
             .names = "{.col}_av"))
# one row per time, containing patch totals and mean individual sizes

The self-thinning plot places mean stem diameter on the horizontal axis and stem density on the vertical axis. Both axes are logarithmic so an approximate power-law relationship appears as a straight line:

patch_total %>%
    ggplot(aes(diameter_stem_av, density)) +
             geom_line() +
             scale_x_log10() +
             scale_y_log10() +
             xlab("stem diameter (m)") +
             ylab("density " ~(m^-2)) +
             ggtitle("Self-thinning plot") +
             theme_classic(base_size = 18)

Leaf area and above-ground biomass provide complementary views of stand development:

patch_total %>%
      ggplot(aes(time, area_leaf)) +
               geom_line() +
               xlab("time (yr)") +
               ylab("leaf area " ~(m^2)) +
               ggtitle("Leaf Area over time") +
               theme_classic(base_size = 18)

patch_total %>%
    ggplot(aes(time, mass_above_ground)) +
             geom_line() +
             xlab("time (yr)") +
             ylab("above-ground biomass (kg)") +
             ggtitle("Above-ground Biomass over time") +
             theme_classic(base_size = 18)

Finally, compare height distributions at selected patch ages. The helper chooses the recorded time nearest each requested age, avoiding assumptions about the number or spacing of solver steps.

target_times <- c(0.5, 2, 5, 9, 14, 28, 56, 98)
recorded_times <- sort(unique(tidy_stand$species$time))
times <- vapply(target_times, function(x) {
  recorded_times[which.min(abs(recorded_times - x))]
}, numeric(1))

# find closest times to a given vector of times for labelling selected size distributions
find_closest <- function(x, at){
  out <- rep("", length(x))

  for(a in at){
    dist <- abs(x-a)
    i <- which(dist == min(dist))
    out[i] <- a
  }
  out
}

data_sub <- tidy_stand$species %>%
  mutate(my_label = find_closest(round(time, 2), at = round(times, 2))) %>%
  na.omit() %>%
  filter(time %in% times)

data_sub2 <- data_sub %>% group_by(time) %>% slice(1) %>% ungroup()

tidy_stand$species %>%
  ggplot(aes(height, density, group=time)) +
  geom_line() +
  geom_line(data = data_sub, col="blue") +
  geom_text(data = data_sub2, aes(label = my_label), nudge_y = 0.05, col="blue") +
  xlab("height (m)") +
  ylab("density " ~(m^-2)) +
  ggtitle("Size distributions over time") +
  theme_classic(base_size = 18) +
  theme(plot.title = element_text(hjust = 0.5))

Young patches contain many similarly sized saplings. As competition becomes asymmetric, density falls and the range of heights broadens because taller plants retain more light. By the oldest highlighted age, few plants remain and their height distribution is narrower near the top of the canopy.