```
using VPL
using Sky
using Plots
using Distributions
import Random
Random.seed!(123456789)
using Base.Threads: @threads
```

# RUE-driven forest

In this example we extend the forest growth model to include PAR interception and a radiation use efficiency to compute the daily growth rate.

The initial setup is as follows:

## Model definition

### Node types

The data types needed to simulate the trees are given in the following module. The difference with respec to the previous model is that Internodes and Leaves have optical properties needed for ray tracing (they are defined as Lambertian surfaces).

```
# Data types
module TreeTypes
using VPL
using Distributions
# Meristem
Base.@kwdef mutable struct Meristem <: VPL.Node
::Int64 = 0 # Age of the meristem
ageend
# Bud
struct Bud <: VPL.Node end
# Node
struct Node <: VPL.Node end
# BudNode
struct BudNode <: VPL.Node end
# Internode (needs to be mutable to allow for changes over time)
Base.@kwdef mutable struct Internode <: VPL.Node
::Int64 = 0 # Age of the internode
age::Float64 = 0.0 # Initial biomass
biomass::Float64 = 0.0 # Internodes
length::Float64 = 0.0 # Internodes
width::Exponential{Float64} = Exponential(5)
sink::Lambertian{1} = Lambertian(τ = 0.1, ρ = 0.05) # Leaf material
materialend
# Leaf
Base.@kwdef mutable struct Leaf <: VPL.Node
::Int64 = 0 # Age of the leaf
age::Float64 = 0.0 # Initial biomass
biomass::Float64 = 0.0 # Leaves
length::Float64 = 0.0 # Leaves
width::Beta{Float64} = Beta(2,5)
sink::Lambertian{1} = Lambertian(τ = 0.1, ρ = 0.05) # Leaf material
materialend
# Graph-level variables -> mutable because we need to modify them during growth
Base.@kwdef mutable struct treeparams
# Variables
::Float64 = 0.0 # Total PAR absorbed by the leaves on the tree (MJ)
PAR::Float64 = 2e-3 # Current total biomass (g)
biomass# Parameters
::Float64 = 5.0 # Radiation use efficiency (g/MJ) -> unrealistic to speed up sim
RUE::Float64 = 1e-3 # Initial biomass of an internode (g)
IB0::Float64 = 0.1e6 # Specific internode weight (g/m3)
SIW::Float64 = 15.0 # Internode shape parameter (length/width)
IS::Float64 = 1e-3 # Initial biomass of a leaf
LB0::Float64 = 100.0 # Specific leaf weight (g/m2)
SLW::Float64 = 3.0 # Leaf shape parameter (length/width)
LS::Float64 = 1/0.5 # Bud break probability coefficient (in 1/m)
budbreak::Int64 = 5 # Number of days between phytomer production
plastochron::Float64 = 15.0 # Number of days that a leaf expands
leaf_expansion::Float64 = 140.0
phyllotaxis::Float64 = 30.0
leaf_angle::Float64 = 45.0
branch_angleend
end
import .TreeTypes
```

### Geometry

The methods for creating the geometry and color of the tree are the same as in the previous example but include the materials for the ray tracer.

```
# Create geometry + color for the internodes
function VPL.feed!(turtle::Turtle, i::TreeTypes.Internode, vars)
# Rotate turtle around the head to implement elliptical phyllotaxis
rh!(turtle, vars.phyllotaxis)
HollowCylinder!(turtle, length = i.length, height = i.width, width = i.width,
= true, color = RGB(0.5,0.4,0.0), material = i.material)
move return nothing
end
# Create geometry + color for the leaves
function VPL.feed!(turtle::Turtle, l::TreeTypes.Leaf, vars)
# Rotate turtle around the arm for insertion angle
ra!(turtle, -vars.leaf_angle)
# Generate the leaf
Ellipse!(turtle, length = l.length, width = l.width, move = false,
= RGB(0.2,0.6,0.2), material = l.material)
color # Rotate turtle back to original direction
ra!(turtle, vars.leaf_angle)
return nothing
end
# Insertion angle for the bud nodes
function VPL.feed!(turtle::Turtle, b::TreeTypes.BudNode, vars)
# Rotate turtle around the arm for insertion angle
ra!(turtle, -vars.branch_angle)
end
```

### Development

The meristem rule is now parameterized by the initial states of the leaves and internodes and will only be triggered every X days where X is the plastochron.

```
# Create right side of the growth rule (parameterized by the initial states
# of the leaves and internodes)
function create_meristem_rule(vleaf, vint)
= Rule(TreeTypes.Meristem,
meristem_rule = mer -> mod(data(mer).age, vars(mer).plastochron) == 0,
lhs = mer -> TreeTypes.Node() +
rhs Bud(),
(TreeTypes.Leaf(biomass = vleaf.biomass,
TreeTypes.= vleaf.length,
length = vleaf.width)) +
width Internode(biomass = vint.biomass,
TreeTypes.= vint.length,
length = vint.width) +
width Meristem())
TreeTypes.end
```

The bud break probability is now a function of distance to the apical meristem rather than the number of internodes. An adhoc traversal is used to compute this length of the main branch a bud belongs to (ignoring the lateral branches).

```
# Compute the probability that a bud breaks as function of distance to the meristem
function prob_break(bud)
# We move to parent node in the branch where the bud was created
= parent(bud)
node # Extract the first internode
= filter(x -> data(x) isa TreeTypes.Internode, children(node))[1]
child = data(child)
data_child # We measure the length of the branch until we find the meristem
= 0.0
distance while !isa(data_child, TreeTypes.Meristem)
# If we encounter an internode, store the length and move to the next node
if data_child isa TreeTypes.Internode
+= data_child.length
distance = children(child)[1]
child = data(child)
data_child # If we encounter a node, extract the next internode
elseif data_child isa TreeTypes.Node
= filter(x -> data(x) isa TreeTypes.Internode, children(child))[1]
child = data(child)
data_child else
error("Should be Internode, Node or Meristem")
end
end
# Compute the probability of bud break as function of distance and
# make stochastic decision
= min(1.0, distance*vars(bud).budbreak)
prob return rand() < prob
end
# Branch rule parameterized by initial states of internodes
function create_branch_rule(vint)
= Rule(TreeTypes.Bud,
branch_rule = prob_break,
lhs = bud -> TreeTypes.BudNode() +
rhs Internode(biomass = vint.biomass,
TreeTypes.= vint.length,
length = vint.width) +
width Meristem())
TreeTypes.end
```

### Light interception

As growth is now dependent on intercepted PAR via RUE, we now need to simulate light interception by the trees. We will use a ray-tracing approach to do so. The first step is to create a scene with the trees and the light sources. As for rendering, the scene can be created from the `forest`

object by simply calling `Scene(forest)`

that will generate the 3D meshes and connect them to their optical properties.

However, we also want to add the soil surface as this will affect the light distribution within the scene due to reflection from the soil surface. This is similar to the customized scene that we created before for rendering, but now for the light simulation.

```
function create_soil()
= Rectangle(length = 21.0, width = 21.0)
soil rotatey!(soil, π/2) # To put it in the XY plane
translate!(soil, Vec(0.0, 10.5, 0.0)) # Corner at (0,0,0)
VPL.return soil
end
function create_scene(forest)
# These are the trees
= Scene(vec(forest))
scene # Add a soil surface
= create_soil()
soil = Lambertian(τ = 0.0, ρ = 0.21)
soil_material add!(scene, mesh = soil, material = soil_material)
# Return the scene
return scene
end
```

Given the scene, we can create the light sources that can approximate the solar irradiance on a given day, location and time of the day using the functions from the Sky package (see package documentation for details). Given the latitude, day of year and fraction of the day (`f = 0`

being sunrise and `f = 1`

being sunset), the function `clear_sky()`

computes the direct and diffuse solar radiation assuming a clear sky. These values may be converted to different wavebands and units using `waveband_conversion()`

. Finally, the collection of light sources approximating the solar irradiance distribution over the sky hemisphere is constructed with the function `sky()`

(this last step requires the 3D scene as input in order to place the light sources adequately).

```
function create_sky(;scene, lat = 52.0*π/180.0, DOY = 182)
# Fraction of the day and day length
= collect(0.1:0.1:0.9)
fs = declination(DOY)
dec = day_length(lat, dec)*3600
DL # Compute solar irradiance
= [clear_sky(lat = lat, DOY = DOY, f = f) for f in fs] # W/m2
temp = getindex.(temp, 1)
Ig = getindex.(temp, 2)
Idir = getindex.(temp, 3)
Idif # Conversion factors to PAR for direct and diffuse irradiance
= waveband_conversion(Itype = :direct, waveband = :PAR, mode = :power)
f_dir = waveband_conversion(Itype = :diffuse, waveband = :PAR, mode = :power)
f_dif # Actual irradiance per waveband
= f_dir.*Idir
Idir_PAR = f_dif.*Idif
Idif_PAR # Create the dome of diffuse light
= sky(scene,
dome = 0.0, # No direct solar radiation
Idir = sum(Idir_PAR)/10*DL, # Daily Diffuse solar radiation
Idif = 1_000_000, # Total number of rays for diffuse solar radiation
nrays_dif = StandardSky, # Angular distribution of solar radiation
sky_model = equal_solid_angles, # Discretization of the sky dome
dome_method = 9, # Number of discretization steps in the zenith angle
ntheta = 12) # Number of discretization steps in the azimuth angle
nphi # Add direct sources for different times of the day
for I in Idir_PAR
push!(dome, sky(scene, Idir = I/10*DL, nrays_dir = 100_000, Idif = 0.0)[1])
end
return dome
end
```

The 3D scene and the light sources are then combined into a `RayTracer`

object, together with general settings for the ray tracing simulation chosen via `RTSettings()`

. The most important settings refer to the Russian roulette system and the grid cloner (see section on Ray Tracing for details). The settings for the Russian roulette system include the number of times a ray will be traced deterministically (`maxiter`

) and the probability that a ray that exceeds `maxiter`

is terminated (`pkill`

). The grid cloner is used to approximate an infinite canopy by replicating the scene in the different directions (`nx`

and `ny`

being the number of replicates in each direction along the x and y axes, respectively). It is also possible to turn on parallelization of the ray tracing simulation by setting `parallel = true`

(currently this uses Julia’s builtin multithreading capabilities).

In addition `RTSettings()`

, an acceleration structure and a splitting rule can be defined when creating the `RayTracer`

object (see ray tracing documentation for details). The acceleration structure allows speeding up the ray tracing by avoiding testing all rays against all objects in the scene.

```
function create_raytracer(scene, sources)
= RTSettings(pkill = 0.9, maxiter = 4, nx = 5, ny = 5, dx = 20.0,
settings = 20.0, parallel = true)
dy RayTracer(scene, sources, settings = settings, acceleration = BVH,
= SAH{3}(5, 10));
rule end
```

The actual ray tracing simulation is performed by calling the `trace!()`

method on the ray tracing object. This will trace all rays from all light sources and update the radiant power absorbed by the different surfaces in the scene inside the `Material`

objects (see `feed!()`

above):

```
function run_raytracer!(forest; DOY = 182)
= create_scene(forest)
scene = create_sky(scene = scene, DOY = DOY)
sources = create_raytracer(scene, sources)
rtobj trace!(rtobj)
return nothing
end
```

The total PAR absorbed for each tree is calculated from the material objects of the different internodes (using `power()`

on the `Material`

object). Note that the `power()`

function returns three different values, one for each waveband, but they are added together as RUE is defined for total PAR.

```
# Run the ray tracer, calculate PAR absorbed per tree and add it to the daily
# total using general weighted quadrature formula
function calculate_PAR!(forest; DOY = 182)
# Reset PAR absorbed by the tree (at the start of a new day)
reset_PAR!(forest)
# Run the ray tracer to compute daily PAR absorption
run_raytracer!(forest, DOY = DOY)
# Add up PAR absorbed by each leaf within each tree
@threads for tree in forest
for l in get_leaves(tree)
vars(tree).PAR += power(l.material)[1]
end
end
return nothing
end
# Reset PAR absorbed by the tree (at the start of a new day)
function reset_PAR!(forest)
for tree in forest
vars(tree).PAR = 0.0
end
return nothing
end
```

### Growth

We need some functions to compute the length and width of a leaf or internode from its biomass

```
function leaf_dims(biomass, vars)
= biomass
leaf_biomass = biomass/vars.SLW
leaf_area = sqrt(leaf_area*4*vars.LS/pi)
leaf_length = leaf_length/vars.LS
leaf_width return leaf_length, leaf_width
end
function int_dims(biomass, vars)
= biomass
int_biomass = biomass/vars.SIW
int_volume = cbrt(int_volume*4*vars.IS^2/pi)
int_length = int_length/vars.IS
int_width return int_length, int_width
end
```

Each day, the total biomass of the tree is updated using a simple RUE formula and the increment of biomass is distributed across the organs proportionally to their relative sink strength (of leaves or internodes).

The sink strength of leaves is modelled with a beta distribution scaled to the `leaf_expansion`

argument that determines the duration of leaf growth, whereas for the internodes it follows a negative exponential distribution. The `pdf`

function computes the probability density of each distribution which is taken as proportional to the sink strength (the model is actually source-limited since we imposed a particular growth rate).

```
sink_strength(leaf, vars) = leaf.age > vars.leaf_expansion ? 0.0 :
pdf(leaf.sink, leaf.age/vars.leaf_expansion)/100.0
plot(0:1:50, x -> sink_strength(TreeTypes.Leaf(age = x), TreeTypes.treeparams()),
= "Age", ylabel = "Sink strength", label = "Leaf") xlabel
```

```
sink_strength(int) = pdf(int.sink, int.age)
plot!(0:1:50, x -> sink_strength(TreeTypes.Internode(age = x)), label = "Internode")
```

Now we need a function that updates the biomass of the tree, allocates it to the different organs and updates the dimensions of said organs. For simplicity, we create the functions `leaves()`

and `internodes()`

that will apply the queries to the tree required to extract said nodes:

```
get_leaves(tree) = apply(tree, Query(TreeTypes.Leaf))
get_internodes(tree) = apply(tree, Query(TreeTypes.Internode))
```

The age of the different organs is updated every time step:

```
function age!(all_leaves, all_internodes, all_meristems)
for leaf in all_leaves
+= 1
leaf.age end
for int in all_internodes
+= 1
int.age end
for mer in all_meristems
+= 1
mer.age end
return nothing
end
```

The daily growth is allocated to different organs proportional to their sink strength.

```
function grow!(tree, all_leaves, all_internodes)
# Compute total biomass increment
= vars(tree)
tvars = max(0.5, tvars.RUE*tvars.PAR/1e6) # Trick to emulate reserves in seedling
ΔB += ΔB
tvars.biomass # Total sink strength
= 0.0
total_sink for leaf in all_leaves
+= sink_strength(leaf, tvars)
total_sink end
for int in all_internodes
+= sink_strength(int)
total_sink end
# Allocate biomass to leaves and internodes
for leaf in all_leaves
+= ΔB*sink_strength(leaf, tvars)/total_sink
leaf.biomass end
for int in all_internodes
+= ΔB*sink_strength(int)/total_sink
int.biomass end
return nothing
end
```

Finally, we need to update the dimensions of the organs. The leaf dimensions are

```
function size_leaves!(all_leaves, tvars)
for leaf in all_leaves
= leaf_dims(leaf.biomass, tvars)
leaf.length, leaf.width end
return nothing
end
function size_internodes!(all_internodes, tvars)
for int in all_internodes
= int_dims(int.biomass, tvars)
int.length, int.width end
return nothing
end
```

### Daily step

All the growth and developmental functions are combined together into a daily step function that updates the forest by iterating over the different trees in parallel.

```
get_meristems(tree) = apply(tree, Query(TreeTypes.Meristem))
function daily_step!(forest, DOY)
# Compute PAR absorbed by each tree
calculate_PAR!(forest, DOY = DOY)
# Grow the trees
@threads for tree in forest
# Retrieve all the relevant organs
= get_leaves(tree)
all_leaves = get_internodes(tree)
all_internodes = get_meristems(tree)
all_meristems # Update the age of the organs
age!(all_leaves, all_internodes, all_meristems)
# Grow the tree
grow!(tree, all_leaves, all_internodes)
= vars(tree)
tvars size_leaves!(all_leaves, tvars)
size_internodes!(all_internodes, tvars)
# Developmental rules
rewrite!(tree)
end
end
```

### Initialization

The trees are initialized on a regular grid with random values for the initial orientation and RUE:

```
= rand(Normal(1.5,0.2), 10, 10)
RUEs histogram(vec(RUEs))
```

```
= [rand()*360.0 for i = 1:2.0:20.0, j = 1:2.0:20.0]
orientations histogram(vec(orientations))
```

`= [Vec(i,j,0) for i = 1:2.0:20.0, j = 1:2.0:20.0]; origins `

The following initalizes a tree based on the origin, orientation and RUE:

```
function create_tree(origin, orientation, RUE)
# Initial state and parameters of the tree
= TreeTypes.treeparams(RUE = RUE)
vars # Initial states of the leaves
= leaf_dims(vars.LB0, vars)
leaf_length, leaf_width = (biomass = vars.LB0, length = leaf_length, width = leaf_width)
vleaf # Initial states of the internodes
= int_dims(vars.LB0, vars)
int_length, int_width = (biomass = vars.IB0, length = int_length, width = int_width)
vint # Growth rules
= create_meristem_rule(vleaf, vint)
meristem_rule = create_branch_rule(vint)
branch_rule = T(origin) + RH(orientation) +
axiom Internode(biomass = vint.biomass,
TreeTypes.= vint.length,
length = vint.width) +
width Meristem()
TreeTypes.= Graph(axiom = axiom, rules = (meristem_rule, branch_rule),
tree = vars)
vars return tree
end
```

## Visualization

As in the previous example, it makes sense to visualize the forest with a soil tile beneath it. Unlike in the previous example, we will construct the soil tile using a dedicated graph and generate a `Scene`

object which can later be merged with the rest of scene generated in daily step:

```
Base.@kwdef struct Soil <: VPL.Node
::Float64
length::Float64
widthend
function VPL.feed!(turtle::Turtle, s::Soil, vars)
Rectangle!(turtle, length = s.length, width = s.width, color = RGB(255/255, 236/255, 179/255))
end
= RA(-90.0) + T(Vec(0.0, 10.0, 0.0)) + # Moves into position
soil_graph Soil(length = 20.0, width = 20.0) # Draws the soil tile
= Scene(Graph(axiom = soil_graph));
soil render(soil, axes = false)
```

And the following function renders the entire scene (notice that we need to use `display()`

to force the rendering of the scene when called within a loop or a function):

```
function render_forest(forest, soil)
= Scene(vec(forest)) # create scene from forest
scene = Scene([scene, soil]) # merges the two scenes
scene display(render(scene))
end
```

## Simulation

We can now create a forest of trees on a regular grid:

```
= create_tree.(origins, orientations, RUEs);
forest render_forest(forest, soil)
= 180
start for i in 1:50
println("Day $i")
daily_step!(forest, i + start)
if mod(i, 5) == 0
render_forest(forest, soil)
end
end
```