Emergent properties

example

An emergent property belongs to the whole patch rather than to any one plant. It arises from many individuals growing, recruiting, dying, and changing one another’s light environment. This example moves from node-density profiles to leaf area and growth, then averages patch profiles across a metapopulation.

Note

A few steps below reach into plant’s internal helpers (plant:::trapezium, plant:::pad_list_to_array). These are not part of the public API and may change or disappear between releases. They are kept here exactly as written because the example depends on them; treat them as implementation detail rather than something to build on.

We’ll work with a single species and an illustrative constant birth rate:

library(plant)

p0 <- scm_base_parameters("FF16")
p <- add_strategies(p0, trait_matrix(0.0825, "lma"), birth_rate = 17.31)
patch <- run_scm(p, refine_schedule = TRUE)$parameters
result <- run_scm(patch, collect = TRUE)

Node-density profiles

Patch profiles change continuously, so plotting every recorded age at equal emphasis is difficult to read. We retain all trajectories in grey and highlight five representative ages:

closest <- function(t, time) {
  which.min(abs(time - t))
}
last <- function(x) {
  x[[length(x)]]
}

times <- c(5, 10, 20, 40, dplyr::last(result$steps$time))
i <- vapply(times, closest, integer(1), result$steps$time)
blues <- c("#DEEBF7", "#C6DBEF", "#9ECAE1", "#6BAED6",
           "#4292C6", "#2171B5", "#08519C", "#08306B")
cols <- colorRampPalette(blues[-(1:2)])(length(i))

species_1   <- dplyr::filter(result$species, species == "1")
height      <- sdd_matrix(species_1, "height")
log_density <- sdd_matrix(species_1, "log_density")

The helper arranges both height and density so rows are nodes and columns are recorded times. Each column can therefore be read as one patch-age profile:

density <- exp(log_density)
matplot(height, density, type="l", lty=1,
        col=util_colour_set_opacity("black", 0.15),
        xlab="Height (m)", ylab="Density (1 / m / m2)", las=1,
        log="y")

Some nodes reach extremely small mathematical densities while their representative individuals continue to grow. This is a continuous-density model, so a value below one does not mean that a literal fraction of a plant is present in a particular physical patch. The highlighted ages make the consequences easier to see.

xlim <- c(0, max(height, na.rm=TRUE) * 1.05)
matplot(height, density, type="l", lty=1,
        col=util_colour_set_opacity("black", 0.15),
        xlab="Height (m)", ylab="Density (1 / m / m2)", las=1,
        log="y", xlim=xlim)
matlines(height[, i], density[, i], col=cols, lty=1, type="l")
points(height[1, i], density[1, i], pch=19, col=cols)
text(height[1, i] + strwidth("x"), density[1, i],
     paste0(round(times), c(" years", rep("", length(times) - 1))),
     adj=c(0, 0))

At five years, self-thinning creates a dip in density. As the stand develops, that dip becomes deeper, broader, and shifts towards taller plants because the low-density nodes still follow the model’s characteristics.

A later wave of recruitment creates another density peak near 4 m at 20 years; by 40 years that cohort has reached roughly 13 m. In old patches, the overall size-density profile changes more slowly, although a narrow gap remains just below the canopy top.

Leaf-area profiles

Which heights contribute most of the patch’s leaf area? The collected output records individual leaf area as competition_effect. Multiplying it by node density gives leaf area density. The canopy model later applies the light extinction coefficient k_I when computing competition; see Canopy and light environments.

species_1 <- species_1 |>
  dplyr::mutate(leaf_area = density * competition_effect)
leaf_area <- sdd_matrix(species_1, "leaf_area")

matplot(height, leaf_area, type="l", lty=1, col="lightgrey",
        xlim=xlim, xlab="Height (m)",
        ylab="Leaf area density (m2 / m2 / m)", las=1)
matlines(height[, i], leaf_area[, i], col=cols, lty=1, type="l")
points(height[1, i], leaf_area[1, i], pch=19, col=cols)
text(height[1, i] + strwidth("x"), leaf_area[1, i],
     paste0(round(times), c(" years", rep("", length(times) - 1))),
     adj=c(0, 0))

Height-growth profiles

Where in the size distribution is height growth concentrated? Unlike the fixed environment in Your first individual, this patch has a shared light environment created by all its plants.

The collected output records states rather than instantaneous rates. We therefore approximate height growth from the change in each node’s height between successive recorded times:

time <- result$steps$time
growth_rate <- t(apply(height, 1, function(hh) c(NA, diff(hh) / diff(time))))

matplot(height, growth_rate, type="l", lty=1, col="lightgrey",
        xlim=xlim, xlab="Height (m)",
        ylab="Height growth rate (m / year)", las=1)
matlines(height[, i], growth_rate[, i], col=cols, lty=1, type="l")
points(height[1, i], growth_rate[1, i], pch=19, col=cols)
text(height[1, i] + strwidth("x"), growth_rate[1, i],
     paste0(round(times), c(" years", rep("", length(times) - 1))),
     adj=c(0, 0))

Average across a metapopulation

The earlier curves each describe a patch of one age. A metapopulation contains patches of many ages, so its average profile weights every patch age by its abundance. Before averaging, we must interpolate each patch profile onto the same set of heights. See the size-structured PDE for the corresponding integral.

hh <- seq_log_range(range(height, na.rm=TRUE), 500)

The next block is reusable numerical scaffolding; readers interested mainly in the ecological result can skip to the resulting average. It interpolates on a log-log scale and returns zero outside the observed height range.

##' Spline interpolation in log-x space
##' @title Spline interpolation in log-x space
##' @param x,y Vectors giving coordinates of points to be
##' interpolated.  The x points should be naturally on a log scale,
##' and for \code{splinefun_loglog} both x and y should be on a log
##' scale.
##' @param ... Additional parameters passed to
##' @rdname splinefun_log
splinefun_loglog <- function(x, y, ...) {
        f <- splinefun(log(x), log(y), ...)
        function(x) {
                exp(f(log(x)))
        }
}

##' Clamp a function to a fixed value outside a domain
##' @param f A function that takes \code{x} as a first argument.
##' @param r Range of values (vector of length 2)
##' @param value (Single) value to use when out of domain.
##' @return A new function
##' @export
clamp_domain <- function(f, r, value = NA_real_) {
        f <- match.fun(f)
        if (length(r) != 2L) {
                stop("Expected length two range")
        }
        if (any(is.na(r)) || r[[2]] < r[[1]]) {
                stop("Values for range must be finite and not decreasing")
        }
        if (length(value) != 1L) {
                stop("value must be length 1")
        }
        function(x, ...) {
                ret <- rep_len(value, length(x))
                i <- x >= r[[1]] & x <= r[[2]]
                if (any(i)) {
                        ret[i] <- f(x[i])
                }
                ret
        }
}

f <- function(height, density, hout) {
  r <- range(height, na.rm=TRUE)
  clamp_domain(splinefun_loglog(height, density), r, 0)(hout)
}

Interpolate every patch-age profile onto the shared height grid:

xx <- lapply(seq_along(result$steps$time),
             function(i) f(height[, i], density[, i], hh))
n_hh <- plant:::pad_list_to_array(xx)

At each height, integrate across patch ages while weighting by patch abundance:

trapezium <- plant:::trapezium
n_av <- apply(n_hh, 1,
              function(x) trapezium(result$steps$time, x * result$steps$patch_density))

The red line is the resulting metapopulation-average density profile:

xlim <- c(0, max(height, na.rm=TRUE) * 1.05)
matplot(height, density, type="l", lty=1,
        col=util_colour_set_opacity("black", 0.15),
        xlab="Height (m)", ylab="Density (1 / m / m2)", las=1,
        log="y", xlim=xlim)
matlines(height[, i], density[, i], col=cols, lty=1, type="l")
points(height[1, i], density[1, i], pch=19, col=cols)
text(height[1, i] + strwidth("x"), density[1, i],
     paste0(round(times), c(" years", rep("", length(times) - 1))),
     adj=c(0, 0))
points(hh, n_av, col="red", type='l')