Light capture in the canopy

Run this notebook

Warning

This area of pyrealm is in active development and this notebook currently contains notes and initial demonstration code.

import numpy as np


from pyrealm.demography.flora import PlantFunctionalType, Flora
from pyrealm.demography.community import Cohorts, Community

from pyrealm.demography.canopy import Canopy, CohortCanopyData

import matplotlib.pyplot as plt

Light capture in a simple model

The big leaf model provides a simple model of the light absorption by a tree. The key variables in determining the light absorption are:

  • The tree has a crown area \(A_c\).

  • The leaf area index \(L\) for the plant functional type is a simple factor that sets the leaf density of the crown, allowing stems with identical crown area to vary in the density of actual leaf surface for light capture. Values of \(L\) are always expressed as the area of leaf surface per square meter of canopy area.

  • The extinction coefficient \(k\) for a PFT sets how much light is absorbed when passing through the leaf surface of that PFT.

The values of \(k\) and \(L\) set the fraction of light absorbed \(f_{abs}\) through the single crown layer following the Beer-Lambert law:

\[ f_{abs} = 1 - e^{-kL} \]

In the big leaf model of a tree canopy, the entire crown area \(A_c\) is exposed to the downwelling photosynthetic photon flux density (PPFD) at the top of the canopy and hence the absorbed radiation (\(I_{abs}\)) is:

\[ I_{abs} = \textrm{PPFD} \cdot A_c (1 - e^{-kL}) \]

Effectively the crown is treated as a single flat area - hence the “big leaf”.

Light capture in a vertically structured canopy

Instead of the big leaf model, the canopy model uses the vertical distribution of leaf area for individual stems, as described in the crown model, to partition all of crown area of individuals in a community into layers. The light capture model needs to account for light absorption through the canopy layers to calculate how much light is absorbed by different individuals within the different layers.

The canopy model defines a set of canopy closure heights (\(z^*_l\) for layers \(l = 1, ..., m\)) that define the canopy layers, but what are these layers exactly? Each layer:

  • has a depth that extends from the closure height of the layer above (or the top of the canopy for the first layer), down to the layer’s closure height.

  • contains projected crown and leaf area from at least one individual in the community.

  • absorbs light, following the projected leaf area of the individuals in the layer and the leaf area index (LAI, \(L\)).

Each successive layer therefore has access to less and less light as it is filtered out by the layers above.

Projected leaf area, leaf area index and layer depth

The projected leaf area and leaf area index sound similar but represent very different things:

  • The projected leaf area is a simple two dimensional measure of how much of the total crown area is found at different heights. Like the projected crown area, the projected leaf area is the cumulative area from the top of the crown down to a height \(z\).

  • To calculate the amount of crown or leaf area between two heights - such as the amount of area within a closed canopy layer - you need to find the difference in projected crown or leaf area between the heights of the layer top and bottom.

  • The total area actual leaf surface within the layer is the leaf area within the layer multiplied by the LAI. Two individuals can have the same crown area and vertical distribution of projected leaf area, but still differ in their LAI. A stem with a larger LAI will have denser canopy that can capture more light.

  • The LAI is effectively the “depth” of leaf area within a layer, because as \(L\) increases, more light is captured within those leaves because \(L\) is included in the Beer-Lambert law.

  • The actual vertical depth of a layer - the distance between two successive canopy closure heights - is not actually used in the light capture model.

For a single individual within the community, the light absorbed in a particular layer \(l\) of the canopy is described as:

\[ I_{abs} = \textrm{PPFD} \cdot T_l \cdot A_{fl} \cdot (1 - e^{-kL}) \]

The first two terms are the PPFD at the top of the canopy and a term \(T_l\) that captures the fraction of the canopy top flux that is transmitted down to layer \(l\).

The last two terms then describe the ability of the individual to capture that downwelling light flux. The Beer-Lambert term is assumed not to change with height - the fraction of light absorbed per metre squared of leaf area does not vary through the canopy. The other term \(A_{fl}\) is the area of foliage of that individual in the layer and is the difference between the projected leaf area at the top and bottom of the layers. That is, given the function from the crown model for the projected leaf area at height \(z\) (\(\tilde{A}_{cp}(z)\)) and a sequence of layer heights \(z_l\):

\[\begin{split} A_{fl} = \begin{cases} \tilde{A}_{cp}(z_l), & \textrm{if } l=1 \\ \tilde{A}_{cp}(z_l) - \tilde{A}_{cp}(z_{l-1}) & \textrm{otherwise} \end{cases} \end{split}\]

Attention

The calculation of \(A_{fl}\) uses the projected leaf area and not the projected crown area.

  • The projected crown area dictates how each stem occupies space in canopy layers: competition for space in the canopy layers use crown area.

  • When the canopy gap fraction (\(f_g\)) is one, then the projected crown area and projected leaf area are equal. However then \(f_g < 1\), the individual tree crown has holes within their crown and projected leaf area is vertically displaced down into the canopy. Because projected leaf area dictates where the actual leaves are located in the canopy, this is what is used to calculate absorbance in the light capture model.

Calculation of light transmission terms

For a given layer \(l\), the community-wide average fraction of light absorbed by the layer (\(a_l\)) can be calculated as the average absorption of each individual, weighted by the total area occupied by the community \(A\). Given a community of \(N\) stems:

\[ a_{l} = \frac{\sum^{i=1}_{N} {A_{fl}}_{i} \cdot (1 - e^{-k_i L_i})}{A} \]

Attention

Note that the calculation of average transmission requires the intermediate calculation of average absorbance in order to account for unoccupied area in the layer. This can arise in incompletely closed layers at the bottom of the canopy or where a canopy or crown gap fraction is set. Because the absorbance of empty space is zero, the calculation of average absorbance implicitly accounts for unfilled space.

The same is not true for transmission, where unfilled area would need to be included in the summation with a transmittance value of one.

The average transmission through the layer can then be calculated as \(t_l = 1 - a_l\) and the cumulative transmission \(T_l\) can be calculated as the cumulative product of layer transmission values after setting \(t_1 = 1\):

\[\begin{split} T_{l} = \begin{cases} 1 & \textrm{if } l=1 \\ T_{l-1} \cdot t_l & \textrm{otherwise} \end{cases} \end{split}\]

The last value of \(T_l\) is the fraction of canopy top radiation passing through the lowest canopy layer and reaching the ground.

Numerical examples of light capture calculation

The code implements a layer by layer simulation of light capture through a simple canopy and then shows that the array calculation used in pyrealm gives the same results as the simulation.

Input data

The input data sets up the stem leaf area within five canopy layers for three cohorts of trees occupying an 8 m2 plot. The cohorts vary in their leaf area index (\(L\)), extinction coefficient (\(k\)) and number of individuals. The crown area of each stem is arranged into vertical layers, representing a set of closed canopy layers as in the perfect plasticity approximation.

# Vertical distribution of projected leaf area for each cohort.
projected_leaf_area = np.array(
    [[6, 2, 0], [10, 6, 0], [13, 9, 1], [13, 11, 4], [13, 11, 5]]
)

# Calculate the leaf area per layer
layer_leaf_area = np.diff(projected_leaf_area, axis=0, prepend=0)
layer_leaf_area

array([[6, 2, 0],
       [4, 4, 0],
       [3, 3, 1],
       [0, 2, 3],
       [0, 0, 1]])

# Light extinction parameters for cohorts.
pft_lai = np.array([2, 1, 2])
pft_par_ext = np.array([0.5, 0.5, 0.6])

n_individuals = np.array([1, 1, 2])
cell_area = 8

# Set an initial light flux in µmol m-2 s-1
initial_ppfd = 1000  # µmol m2 s1

Just to show the layer structure in more detail, if we multiply the stem leaf area by the number of individuals, we get the leaf area from each cohort in each layer, and then the sum of values within layers shows the leaf area occupying the available cell area within each layer except for the last.

cohort_leaf_area = layer_leaf_area * n_individuals
cohort_leaf_area

array([[6, 2, 0],
       [4, 4, 0],
       [3, 3, 2],
       [0, 2, 6],
       [0, 0, 2]])

cohort_leaf_area.sum(axis=1)

array([8, 8, 8, 8, 2])

Layer-wise simulation

The loop below iterates over the canopy layers from the top layer to the ground. At each layer:

  • The absorption fraction of each layer for each cohort is calculated following the Beer-Lambert law as \(1- e^{-kL}\) and used to calculate:

    • the actual PPFD absorbed per metre squared

    • the fraction of the initial PPFD absorbed in that layer by the leaf area of each cohort.

  • The stem per m2 values are scaled up to the community level using the stem leaf area in each cohort and the number of individuals and used to calculate:

    • the total flux captured across the cell in that layer

    • the average captured flux per m2 and hence

    • the remaining PPFD reaching the layer below and

    • the average fraction of initial PPFD transmitted through the layer

# Set the initial PPFD
ppfd = initial_ppfd

# Create stores for variable created and updated in the loop
simulated_transmission = np.empty((6, 1))
simulated_transmission[0, 0] = 1
simulated_per_stem_f_abs = np.empty((5, 3))
simulated_total_capture = 0

# Loop over the layers
for layer in np.arange(5):

    # Calculate the per stem light capture per m2
    per_stem_light_abs = ppfd * (1 - np.exp(-pft_par_ext * pft_lai))
    simulated_per_stem_f_abs[layer, :] = per_stem_light_abs / initial_ppfd

    # Calculate the total light captured across the layer
    cohort_capture = per_stem_light_abs * layer_leaf_area[layer,] * n_individuals

    # Scale that back down to give the average captured flux per m2
    average_captured_ppfd = cohort_capture.sum() / cell_area
    simulated_total_capture += cohort_capture.sum()

    # Calculate the remaining flux transmitted to the next layer
    ppfd -= average_captured_ppfd
    simulated_transmission[layer + 1, 0] = ppfd / initial_ppfd

The resulting transmission starts with the fraction of light reaching the first layer (always 1) and ends with the fraction of light passing through the last layer to the ground.

simulated_transmission

array([[1.        ],
       [0.42754225],
       [0.20830074],
       [0.09179862],
       [0.03465658],
       [0.02860203]])

The per stem fraction absorbed gives the fraction of the canopy top incident light absorbed by stem leaf area for each cohort in each layer.

Critically these values do not assume any stem leaf area is actually present in the layer. They are what the leaves could achieve if present and so need to be multiplied by stem leaf area to give actual absorbance.

simulated_per_stem_f_abs

array([[0.63212056, 0.39346934, 0.69880579],
       [0.27025824, 0.16822477, 0.298769  ],
       [0.13167118, 0.08195996, 0.14556176],
       [0.0580278 , 0.03611994, 0.06414941],
       [0.02190714, 0.0136363 , 0.02421822]])

Lastly - as a check sum - the total capture shows how much of the incident light across the whole cell (PPFD x cell area) is captured by leaves.

simulated_total_capture

np.float64(7771.183792242919)

Matrix approach

The code below calculates the same quantities more efficiently using array based calculation.

First, we can calculate an array of Beer-Lambert absorption coefficients for the stems in each cohort: this will be constant across array layers.

stem_absorption = 1 - np.exp(-pft_par_ext * pft_lai)
stem_absorption

array([0.63212056, 0.39346934, 0.69880579])

Next, we calculate the average absorption within each layer. We find the sum of the product of the per stem absorptions and the total cohort leaf areas and then divide by the cell area to give the average fraction absorbed (\(a_l\)) for a given square meter.

# Calculate the layer average absorption
average_layer_absorption = (stem_absorption * cohort_leaf_area).sum(axis=1) / cell_area
average_layer_absorption

array([0.57245775, 0.51279495, 0.55929766, 0.62247168, 0.17470145])

We can also calculate the average leaf area index within the layer - this is useful for easy calculation of the total leaf surface within the layer across cohorts as the product of cell area and average leaf area index.

cohort_leaf_area_index = np.broadcast_to(pft_lai, layer_leaf_area.shape)

average_layer_lai = (cohort_leaf_area_index * cohort_leaf_area).sum(axis=1) / cell_area

average_layer_lai

array([1.75 , 1.5  , 1.625, 1.75 , 0.5  ])

We can now calculate average transmission through each layers (\(t_l = 1 - a_l\)):

average_layer_transmission = 1 - average_layer_absorption
average_layer_transmission

array([0.42754225, 0.48720505, 0.44070234, 0.37752832, 0.82529855])

And hence the cumulative product of the average per layer transmissions to give the fraction of downwelling light flux transmitted to each layer (\(T_l\)):

transmission_profile = np.cumprod(np.concat([[1], average_layer_transmission]))[:, None]
transmission_profile

array([[1.        ],
       [0.42754225],
       [0.20830074],
       [0.09179862],
       [0.03465658],
       [0.02860203]])

We can then check that the simulated and matrix based calculations are the same:

if np.allclose(simulated_transmission, transmission_profile):
    print("\U00002705 Transmission matches")

✅ Transmission matches

We can also calculate the average fraction of absorbed photosynthetic radiation (\(f_{APAR}\)) absorbed by each layer. This is simply the differences between the transmission profiles

fapar = np.diff(-transmission_profile, axis=0)
fapar

array([[0.57245775],
       [0.2192415 ],
       [0.11650212],
       [0.05714204],
       [0.00605455]])

Note that the sum of the average layer \(f_{APAR}\) and the last value in the transmission profile, which is the light fraction transmitted to the ground, equals one: all light is absorbed by layers or reaches the ground.

fapar.sum() + transmission_profile[-1]

array([1.])

Next, we can calculate a matrix of per stem \(f_{APAR}\) values as the product of the per layer stem absorptions and the per layer average transmission.

per_stem_f_abs = stem_absorption * transmission_profile[:-1]
per_stem_f_abs

array([[0.63212056, 0.39346934, 0.69880579],
       [0.27025824, 0.16822477, 0.298769  ],
       [0.13167118, 0.08195996, 0.14556176],
       [0.0580278 , 0.03611994, 0.06414941],
       [0.02190714, 0.0136363 , 0.02421822]])

And again check that we get the same solution:

if np.allclose(simulated_per_stem_f_abs, per_stem_f_abs):
    print("\U00002705 Per stem fraction absorbed matches")

✅ Per stem fraction absorbed matches

We can check that matches the average layer \(f_{APAR}\): if we take the per layer sums of the per stem \(f_{APAR}\), weighted by the cohort leaf areas, and divide through by the available cell area, we recreate the average per layer values:

(per_stem_f_abs * cohort_leaf_area).sum(axis=1, keepdims=True) / cell_area

array([[0.57245775],
       [0.2192415 ],
       [0.11650212],
       [0.05714204],
       [0.00605455]])

Last, we can calculate the total absorption in each layer for each stem. This is the key result for use in most cases: how much light flux is absorbed by each stem in the leaf area of each layer.

per_stem_absorption = per_stem_f_abs * layer_leaf_area * initial_ppfd
per_stem_absorption

array([[3792.72335297,  786.93868057,    0.        ],
       [1081.03297337,  672.89906161,    0.        ],
       [ 395.01354324,  245.87986594,  145.56176378],
       [   0.        ,   72.23988832,  192.44823009],
       [   0.        ,    0.        ,   24.21821925]])

Lastly, we can verify that the total absorption across all individuals within the community is the same as the simulation:

total_capture = (per_stem_absorption * n_individuals).sum()
total_capture

np.float64(7771.183792242919)

if np.allclose(simulated_total_capture, total_capture):
    print("\U00002705 Total capture matches")

✅ Total capture matches

Implementation in pyrealm

The calculation above can be directly repeated using the CohortCanopyData class: this data class is not typically used directly but is designed for use with the input data defined above.

canopy_light_model = CohortCanopyData(
    projected_leaf_area=projected_leaf_area,
    lai=pft_lai,
    par_ext=pft_par_ext,
    n_individuals=n_individuals,
    cell_area=cell_area,
)

/home/docs/checkouts/readthedocs.org/user_builds/pyrealm/checkouts/latest/pyrealm/core/experimental.py:72: ExperimentalFeatureWarning: 'Be aware that CohortCanopyData is an experimental feature of pyrealm and the implementation and API may change within major versions.'
  warn(qualname, ExperimentalFeatureWarning)
/home/docs/checkouts/readthedocs.org/user_builds/pyrealm/checkouts/latest/pyrealm/core/experimental.py:72: ExperimentalFeatureWarning: 'Be aware that CommunityCanopyData is an experimental feature of pyrealm and the implementation and API may change within major versions.'
  warn(qualname, ExperimentalFeatureWarning)

The light transmission profile and the final light transmission to the ground match the calculations above.

canopy_light_model.community_data.transmission_profile

array([1.        , 0.42754225, 0.20830074, 0.09179862, 0.03465658])

canopy_light_model.community_data.transmission_to_ground

np.float64(0.028602025969635302)

Similarly, the per stem light absorption matches:

initial_ppfd * canopy_light_model.stem_leaf_area * canopy_light_model.fapar

array([[3792.72335297,  786.93868057,    0.        ],
       [1081.03297337,  672.89906161,    0.        ],
       [ 395.01354324,  245.87986594,  145.56176378],
       [   0.        ,   72.23988832,  192.44823009],
       [   0.        ,    0.        ,   24.21821925]])

Calculation of light capture for a community

The Canopy model automatically calculates the light capture model for a community. The code below uses the final example community from the documentation of the vertical canopy model calculations.

# Define gappy PFTs
gappy_short_pft = PlantFunctionalType(
    name="short",
    h_max=15,
    m=1.5,
    n=1.5,
    f_g=0.1,
    ca_ratio=380,
)
gappy_tall_pft = PlantFunctionalType(
    name="tall", h_max=30, m=3, n=1.5, par_ext=0.6, f_g=0.1, ca_ratio=500
)

# Create the flora
gappy_flora = Flora([gappy_short_pft, gappy_tall_pft])

# Define community with three cohorts
gappy_community = Community(
    flora=gappy_flora,
    cell_area=150,
    cell_id=1,
    cohorts=Cohorts(
        dbh_values=np.array([0.1, 0.20, 0.5]),
        n_individuals=np.array([7, 3, 2]),
        pft_names=np.array(["short", "short", "tall"]),
    ),
)

/home/docs/checkouts/readthedocs.org/user_builds/pyrealm/checkouts/latest/pyrealm/core/experimental.py:72: ExperimentalFeatureWarning: 'Be aware that Cohorts is an experimental feature of pyrealm and the implementation and API may change within major versions.'
  warn(qualname, ExperimentalFeatureWarning)
/home/docs/checkouts/readthedocs.org/user_builds/pyrealm/checkouts/latest/pyrealm/core/experimental.py:72: ExperimentalFeatureWarning: 'Be aware that Community is an experimental feature of pyrealm and the implementation and API may change within major versions.'
  warn(qualname, ExperimentalFeatureWarning)
/home/docs/checkouts/readthedocs.org/user_builds/pyrealm/checkouts/latest/pyrealm/core/experimental.py:72: ExperimentalFeatureWarning: 'Be aware that StemAllometry is an experimental feature of pyrealm and the implementation and API may change within major versions.'
  warn(qualname, ExperimentalFeatureWarning)

We can use the community object directly to calculate the big leaf approximation of light capture.

(
    initial_ppfd
    * gappy_community.stem_allometry.crown_area
    * (
        1
        - np.exp(-gappy_community.stem_traits.par_ext * gappy_community.stem_traits.lai)
    )
).round(2)

array([[ 1233.34,  3605.  , 28683.97]])

If the cell area of the community is set to be greater than the total crown area of the community, a canopy model fitted to the community gives the same estimates: there is a single canopy layer with no shading of any of the stem canopies.

# Calculate the canopy profile across vertical heights
gappy_canopy_ppa = Canopy(community=gappy_community, fit_ppa=True)

# Layer closure heights
gappy_canopy_ppa.heights

/home/docs/checkouts/readthedocs.org/user_builds/pyrealm/checkouts/latest/pyrealm/core/experimental.py:72: ExperimentalFeatureWarning: 'Be aware that Canopy is an experimental feature of pyrealm and the implementation and API may change within major versions.'
  warn(qualname, ExperimentalFeatureWarning)
/home/docs/checkouts/readthedocs.org/user_builds/pyrealm/checkouts/latest/pyrealm/core/experimental.py:72: ExperimentalFeatureWarning: 'Be aware that CrownProfile is an experimental feature of pyrealm and the implementation and API may change within major versions.'
  warn(qualname, ExperimentalFeatureWarning)

array([[0.]])

The canopy model then contains the data required to calculate the absorbed irradiance within each layer for a stem in each cohort.

i_abs = (
    initial_ppfd
    * gappy_canopy_ppa.cohort_data.fapar
    * gappy_canopy_ppa.cohort_data.stem_leaf_area
)
i_abs.sum(axis=0).round(2)

array([ 1233.34,  3605.  , 28683.97])

The example code below places the same community of cohorts into a much smaller area: the tree crowns cannot all fit into a single layer and so are arranged under the PPA model into layers.

gappy_community = Community(
    flora=gappy_flora,
    cell_area=32,
    cell_id=1,
    cohorts=Cohorts(
        dbh_values=np.array([0.1, 0.20, 0.5]),
        n_individuals=np.array([7, 3, 2]),
        pft_names=np.array(["short", "short", "tall"]),
    ),
)

# Recalculate the canopy profile across vertical heights
gappy_canopy_ppa = Canopy(community=gappy_community, fit_ppa=True)

/home/docs/checkouts/readthedocs.org/user_builds/pyrealm/checkouts/latest/pyrealm/demography/crown.py:72: RuntimeWarning: invalid value encountered in power
  q_z = m * n * z_over_height ** (n - 1) * (1 - z_over_height**n) ** (m - 1)

The PPA solution for this community finds four layers and we can show:

  • the canopy closure heights of those layers,

  • the average \(f_{APAR}\) of each layer, and

  • the amount of leaf area present in each layer for each of the three cohorts.

# Layer closure heights
gappy_canopy_ppa.heights

array([[15.02533528],
       [11.59484225],
       [ 8.95730042],
       [ 0.        ]])

gappy_canopy_ppa.community_data.average_layer_fapar

array([0.5943422 , 0.23999871, 0.09442663, 0.0459118 ])

gappy_canopy_ppa.cohort_data.stem_leaf_area.round(2)

array([[ 0.  ,  0.  , 14.4 ],
       [ 0.  ,  0.44, 13.74],
       [ 0.  ,  3.88,  8.58],
       [ 2.08,  1.76,  6.71]])

All of the first layer is occupied by the crowns of the two individuals in the tallest cohort and the rest of the leaf area is shaded, reducing the overall light absorbed by each stem below the big leaf prediction. This is most marked in the smallest stems, which only have leaf area in the very lowest layer.

i_abs = (
    initial_ppfd
    * gappy_canopy_ppa.cohort_data.fapar
    * gappy_canopy_ppa.cohort_data.stem_leaf_area
)
i_abs.sum(axis=0).round(2)

array([   87.85,   561.39, 14445.29])