library(plant)
library(dplyr)
library(tidyr)
library(ggplot2)
library(patchwork)Example analysis
example
This worked example follows a complete analysis from model setup to figures. We will create two TF24 strategies, run them together in one patch, turn the size distributions into patch totals, and interpolate the output to requested times and heights.
NoteCustomising the
TF24 soil-water environment
The TF24 model couples plant growth to a multi-layer soil-water balance. Environment("TF24") already provides five soil layers, an initial moisture state, and constant rainfall. The helper below shows how to replace those defaults with 15 layers and seasonal rainfall.
TF24’s hydraulic calculations make long runs expensive, so this example stops at a patch age of 25 years. Choose a lifetime appropriate to your scientific question for a production analysis.
Set up the analysis
Console logging is optional, but it makes progress visible during the slower TF24 run:
plant_log_console()Define the shortened patch lifetime and a function that constructs the same custom environment every time the model is run:
max_patch_lifetime <- 25
make_tf24_env <- function(max_patch_lifetime, n_layers = 15, theta0 = 0.2) {
env <- Environment("TF24")
env$set_soil_number_of_depths(n_layers)
env$set_soil_water_state(rep(theta0, n_layers))
# seasonal rainfall around a mean of 1.0
x <- seq(0, max_patch_lifetime, length.out = 1000)
y <- 0.25 * sin(2 * pi * x) + 1.0
env$extrinsic_drivers_set_variable("rainfall", x = x, y = y)
env
}Run the smallest version
A minimal run needs base parameters, at least one strategy, and an environment. The first call refines the node-introduction schedule; the second runs the patch and collects its state:
p0 <- scm_base_parameters("TF24", "TF24_Env")
p0$max_patch_lifetime <- max_patch_lifetime
traits <- trait_matrix(c(0.07), c("lma"))
p1 <- add_strategies(p0, traits)
env <- make_tf24_env(max_patch_lifetime)
results <- run_scm(p1, env = env, refine_schedule = TRUE)$parameters
results <- run_scm(results, env = env, collect = TRUE)For repeated analyses, helper functions keep assumptions in one place. The first helper starts from TF24 defaults and changes four strategy parameters:
base_parameters <- function() {
p0 <- scm_base_parameters("TF24", "TF24_Env")
p0$max_patch_lifetime <- max_patch_lifetime
p0$strategy_default$pars$hmat <- 15
p0$strategy_default$pars$rho <- 700
p0$strategy_default$pars$a_l1 <- 2.17
p0$strategy_default$pars$a_l2 <- 0.5
p0
}The run helper below combines three kinds of input:
- biological choices: traits, birth rates, and
B_lf1; - location: latitude, used by the TF24 hyperparameterisation;
- a numerical choice: whether to refine the node-introduction schedule.
run_mypatch <- function(
traits = trait_matrix(c(0.07), c("lma")),
birth_rate = 10,
B_lf1 = 0.009,
p0 = base_parameters(),
optimise_schedule = FALSE,
latitude = 28.182
) {
hyper_par_fn = make_TF24_hyperpar(B_lf1 = B_lf1, latitude = latitude)
birth_rate <- rep_len(birth_rate, nrow(traits))
p1 <- add_strategies(p0, traits, hyper_par_fn, birth_rate = birth_rate)
env <- make_tf24_env(p0$max_patch_lifetime)
if(optimise_schedule)
result <- run_scm(p1, env = env, refine_schedule = TRUE)$parameters
else
result <- p1
# gather outputs at each time step
run_scm(result, env = env, collect = TRUE)
}Run two strategies together
Run a patch for two strategies. Here, we provide two LMA values, representing thin leaves (0.07) and thicker leaves (0.24). To run one species, supply a one-row trait matrix; add_strategies() requires at least one strategy before the patch can be simulated.
results <-
run_mypatch(traits=trait_matrix(c(0.07, 0.24), c("lma")))Because run_mypatch() calls run_scm(..., collect = TRUE), its result is already organised for analysis. The species table has one row per node at each recorded step; steps and env describe patch time and environment.
results_tidy <- results
results_tidy$steps
# A tibble: 102 × 3
step time patch_density
<int> <dbl> <dbl>
1 1 0 0.140
2 2 0.00001 0.140
3 3 0.00002 0.140
4 4 0.00003 0.140
5 5 0.00004 0.140
6 6 0.00005 0.140
7 7 0.00006 0.140
8 8 0.00007 0.140
9 9 0.00008 0.140
10 10 0.0000953 0.140
# ℹ 92 more rows
$n_spp
[1] 2
$species
# A tibble: 10,506 × 25
species time step patch_density node density log_density height mortality
<chr> <dbl> <int> <dbl> <int> <dbl> <dbl> <dbl> <dbl>
1 1 0 1 0.140 1 22.3 3.10 0.0351 0.000699
2 1 1e-5 2 0.140 1 22.3 3.10 0.0351 0.000699
3 1 1e-5 2 0.140 2 22.3 3.10 0.0351 0.000699
4 1 2e-5 3 0.140 1 22.3 3.10 0.0351 0.000699
5 1 2e-5 3 0.140 2 22.3 3.10 0.0351 0.000699
6 1 2e-5 3 0.140 3 22.3 3.10 0.0351 0.000699
7 1 3e-5 4 0.140 1 22.3 3.10 0.0351 0.000699
8 1 3e-5 4 0.140 2 22.3 3.10 0.0351 0.000699
9 1 3e-5 4 0.140 3 22.3 3.10 0.0351 0.000699
10 1 3e-5 4 0.140 4 22.3 3.10 0.0351 0.000699
# ℹ 10,496 more rows
# ℹ 16 more variables: fecundity <dbl>, area_heartwood <dbl>,
# mass_heartwood <dbl>, storage <dbl>,
# offspring_produced_survival_weighted <dbl>, competition_effect <dbl>,
# height_inverse <dbl>, net_mass_production_dt <dbl>, root_mass <dbl>,
# opt_psi_stem <dbl>, opt_root_psi <dbl>, transpiration <dbl>, E_up_ <dbl>,
# profit <dbl>, stom_cond_CO2 <dbl>, assimilation <dbl>
$env
$env$light_availability
# A tibble: 5,482 × 5
time step patch_density height light_availability
<dbl> <int> <dbl> <dbl> <dbl>
1 0 1 0.140 0 1
2 0 1 0.140 0.00110 1
3 0 1 0.140 0.00219 1
4 0 1 0.140 0.00329 1
5 0 1 0.140 0.00438 1
6 0 1 0.140 0.00548 1
7 0 1 0.140 0.00658 1
8 0 1 0.140 0.00767 1
9 0 1 0.140 0.00877 1
10 0 1 0.140 0.00986 1
# ℹ 5,472 more rows
$env$soil_moist
# A tibble: 1,530 × 4
time step patch_density soil_moist
<dbl> <int> <dbl> <dbl>
1 0 1 0.140 0.2
2 0 1 0.140 0.2
3 0 1 0.140 0.2
4 0 1 0.140 0.2
5 0 1 0.140 0.2
6 0 1 0.140 0.2
7 0 1 0.140 0.2
8 0 1 0.140 0.2
9 0 1 0.140 0.2
10 0 1 0.140 0.2
# ℹ 1,520 more rows
$env$soil_depth
# A tibble: 1,530 × 4
time step patch_density soil_depth
<dbl> <int> <dbl> <dbl>
1 0 1 0.140 0.1
2 0 1 0.140 0.2
3 0 1 0.140 0.3
4 0 1 0.140 0.4
5 0 1 0.140 0.5
6 0 1 0.140 0.6
7 0 1 0.140 0.7
8 0 1 0.140 0.8
9 0 1 0.140 0.9
10 0 1 0.140 1
# ℹ 1,520 more rows
$env$soil_moist_cumulative_flux
# A tibble: 102 × 7
time step patch_density sum_rainfall sum_infiltration sum_drainage
<dbl> <int> <dbl> <dbl> <dbl> <dbl>
1 0 1 0.140 0 0 0
2 0.00001 2 0.140 0.0000100 0.00000998 0.00000000758
3 0.00002 3 0.140 0.0000200 0.0000200 0.0000000152
4 0.00003 4 0.140 0.0000300 0.0000299 0.0000000227
5 0.00004 5 0.140 0.0000400 0.0000399 0.0000000303
6 0.00005 6 0.140 0.0000500 0.0000499 0.0000000379
7 0.00006 7 0.140 0.0000600 0.0000599 0.0000000455
8 0.00007 8 0.140 0.0000700 0.0000698 0.0000000530
9 0.00008 9 0.140 0.0000800 0.0000798 0.0000000606
10 0.0000953 10 0.140 0.0000953 0.0000950 0.0000000722
# ℹ 92 more rows
# ℹ 1 more variable: sum_resource_depletion <dbl>
$offspring_production
[1] 8.908685e+00 6.607032e-22
$net_reproduction_ratios
[1] 8.908685e-01 6.607032e-23
$p
$patch_area
[1] 1
$n_patches
[1] 1
$patch_type
[1] "meta-population"
$max_patch_lifetime
[1] 25
$strategies
$strategies[[1]]
$pars
$lma
[1] 0.07
$rho
[1] 700
$hmat
[1] 15
$omega
[1] 3.8e-05
$eta
[1] 12
$theta
[1] 0.0002141786
$a_l1
[1] 2.17
$a_l2
[1] 0.5
$a_r1
[1] 0.07
$a_b1
[1] 0.17
$r_s
[1] 5.731429
$r_b
[1] 11.46286
$r_r
[1] 217
$r_l
[1] 640.9011
$a_y
[1] 0.7
$a_bio
[1] 0.0245
$k_l
[1] 2.699283
$k_b
[1] 0.2
$k_s
[1] 0.2
$k_r
[1] 1
$a_p1
[1] 151.1778
$a_p2
[1] 0.2047162
$a_f3
[1] 0.000114
$a_f1
[1] 1
$a_f2
[1] 50
$S_D
[1] 0.25
$a_d0
[1] 0.1
$d_I
[1] 0.01
$a_dG1
[1] 5.5
$a_dG2
[1] 20
$a_st1
[1] 0.1
$a_st2
[1] 0.1
$a_st3
[1] 0.8
$k_I
[1] 0.5
$vcmax_25
[1] 96
$p_50
[1] 2.888726
$K_s
[1] 1
$c
[1] 2.04
$b
[1] 3.457268
$psi_crit
[1] 5.91988
$beta1
[1] 20000
$beta2
[1] 1.5
$g1_TF24
[1] 7.5
$jmax_25
[1] 157.44
$a
[1] 0.3
$curv_fact_elec_trans
[1] 0.7
$curv_fact_colim
[1] 0.99
$var_sapwood_volume_cost
[1] 1
$nmass_l
[1] 0.02392789
$nmass_s
[1] 0.00198
$nmass_b
[1] 0.0034
$nmass_r
[1] 0.00335
$dmass_dN
[1] 0
$root_depth_shape_eta
[1] 0.2
$recruitment_decay
[1] 0
attr(,"class")
[1] "TF24_Pars"
$control
$function_integration_rule
[1] 21
$shading_model
[1] ""
$ppa_layer_optical_depth
[1] 0.5
$ppa_layer_smoothing
[1] 0.3
$offspring_production_tol
[1] 1e-08
$offspring_production_iterations
[1] 1000
$node_gradient_eps
[1] 1e-06
$node_gradient_direction
[1] -1
$node_gradient_richardson
[1] FALSE
$node_gradient_richardson_depth
[1] 4
$ode_step_size_initial
[1] 1e-06
$ode_step_size_min
[1] 1e-06
$ode_step_size_max
[1] 5
$ode_tol_rel
[1] 1e-04
$ode_tol_abs
[1] 1e-04
$ode_a_y
[1] 1
$ode_a_dydt
[1] 0
$fixed_time_step
[1] 0
$schedule_nsteps
[1] 20
$schedule_eps
[1] 0.02
$schedule_verbose
[1] FALSE
$save_RK45_cache
[1] FALSE
$GSS_tol_abs
[1] 0.001
$vulnerability_curve_ncontrol
[1] 100
$ci_abs_tol
[1] 0.001
$ci_niter
[1] 1000
attr(,"class")
[1] "Control"
$collect_all_auxiliary
[1] FALSE
$birth_rate_x
numeric(0)
$birth_rate_y
[1] 10
$is_variable_birth_rate
[1] FALSE
attr(,"class")
[1] "TF24_Strategy"
$strategies[[2]]
$pars
$lma
[1] 0.24
$rho
[1] 700
$hmat
[1] 15
$omega
[1] 3.8e-05
$eta
[1] 12
$theta
[1] 0.0002141786
$a_l1
[1] 2.17
$a_l2
[1] 0.5
$a_r1
[1] 0.07
$a_b1
[1] 0.17
$r_s
[1] 5.731429
$r_b
[1] 11.46286
$r_r
[1] 217
$r_l
[1] 320.8045
$a_y
[1] 0.7
$a_bio
[1] 0.0245
$k_l
[1] 0.32825
$k_b
[1] 0.2
$k_s
[1] 0.2
$k_r
[1] 1
$a_p1
[1] 151.1778
$a_p2
[1] 0.2047162
$a_f3
[1] 0.000114
$a_f1
[1] 1
$a_f2
[1] 50
$S_D
[1] 0.25
$a_d0
[1] 0.1
$d_I
[1] 0.01
$a_dG1
[1] 5.5
$a_dG2
[1] 20
$a_st1
[1] 0.1
$a_st2
[1] 0.1
$a_st3
[1] 0.8
$k_I
[1] 0.5
$vcmax_25
[1] 96
$p_50
[1] 2.888726
$K_s
[1] 1
$c
[1] 2.04
$b
[1] 3.457268
$psi_crit
[1] 5.91988
$beta1
[1] 20000
$beta2
[1] 1.5
$g1_TF24
[1] 7.5
$jmax_25
[1] 157.44
$a
[1] 0.3
$curv_fact_elec_trans
[1] 0.7
$curv_fact_colim
[1] 0.99
$var_sapwood_volume_cost
[1] 1
$nmass_l
[1] 0.01335397
$nmass_s
[1] 0.00198
$nmass_b
[1] 0.0034
$nmass_r
[1] 0.00335
$dmass_dN
[1] 0
$root_depth_shape_eta
[1] 0.2
$recruitment_decay
[1] 0
attr(,"class")
[1] "TF24_Pars"
$control
$function_integration_rule
[1] 21
$shading_model
[1] ""
$ppa_layer_optical_depth
[1] 0.5
$ppa_layer_smoothing
[1] 0.3
$offspring_production_tol
[1] 1e-08
$offspring_production_iterations
[1] 1000
$node_gradient_eps
[1] 1e-06
$node_gradient_direction
[1] -1
$node_gradient_richardson
[1] FALSE
$node_gradient_richardson_depth
[1] 4
$ode_step_size_initial
[1] 1e-06
$ode_step_size_min
[1] 1e-06
$ode_step_size_max
[1] 5
$ode_tol_rel
[1] 1e-04
$ode_tol_abs
[1] 1e-04
$ode_a_y
[1] 1
$ode_a_dydt
[1] 0
$fixed_time_step
[1] 0
$schedule_nsteps
[1] 20
$schedule_eps
[1] 0.02
$schedule_verbose
[1] FALSE
$save_RK45_cache
[1] FALSE
$GSS_tol_abs
[1] 0.001
$vulnerability_curve_ncontrol
[1] 100
$ci_abs_tol
[1] 0.001
$ci_niter
[1] 1000
attr(,"class")
[1] "Control"
$collect_all_auxiliary
[1] FALSE
$birth_rate_x
numeric(0)
$birth_rate_y
[1] 10
$is_variable_birth_rate
[1] FALSE
attr(,"class")
[1] "TF24_Strategy"
$strategy_default
$pars
$lma
[1] 0.1978791
$rho
[1] 700
$hmat
[1] 15
$omega
[1] 3.8e-05
$eta
[1] 12
$theta
[1] 0.0002141786
$a_l1
[1] 2.17
$a_l2
[1] 0.5
$a_r1
[1] 0.07
$a_b1
[1] 0.17
$r_s
[1] 6.598684
$r_b
[1] 13.19737
$r_r
[1] 217
$r_l
[1] 198.4545
$a_y
[1] 0.7
$a_bio
[1] 0.0245
$k_l
[1] 0.4565855
$k_b
[1] 0.2
$k_s
[1] 0.2
$k_r
[1] 1
$a_p1
[1] 151.1778
$a_p2
[1] 0.2047162
$a_f3
[1] 0.000114
$a_f1
[1] 1
$a_f2
[1] 50
$S_D
[1] 0.25
$a_d0
[1] 0.1
$d_I
[1] 0.01
$a_dG1
[1] 5.5
$a_dG2
[1] 20
$a_st1
[1] 0.1
$a_st2
[1] 0.1
$a_st3
[1] 0.8
$k_I
[1] 0.5
$vcmax_25
[1] 96
$p_50
[1] 1.85
$K_s
[1] 1
$c
[1] 1.089985
$b
[1] 2.589437
$psi_crit
[1] 7.085493
$beta1
[1] 20000
$beta2
[1] 1.5
$g1_TF24
[1] 7.5
$jmax_25
[1] 157.44
$a
[1] 0.3
$curv_fact_elec_trans
[1] 0.7
$curv_fact_colim
[1] 0.99
$var_sapwood_volume_cost
[1] 1
$nmass_l
[1] 0.013
$nmass_s
[1] 0.00198
$nmass_b
[1] 0.0034
$nmass_r
[1] 0.00335
$dmass_dN
[1] 0
$root_depth_shape_eta
[1] 0.2
$recruitment_decay
[1] 0
attr(,"class")
[1] "TF24_Pars"
$control
$function_integration_rule
[1] 21
$shading_model
[1] ""
$ppa_layer_optical_depth
[1] 0.5
$ppa_layer_smoothing
[1] 0.3
$offspring_production_tol
[1] 1e-08
$offspring_production_iterations
[1] 1000
$node_gradient_eps
[1] 1e-06
$node_gradient_direction
[1] -1
$node_gradient_richardson
[1] FALSE
$node_gradient_richardson_depth
[1] 4
$ode_step_size_initial
[1] 1e-06
$ode_step_size_min
[1] 1e-06
$ode_step_size_max
[1] 5
$ode_tol_rel
[1] 1e-04
$ode_tol_abs
[1] 1e-04
$ode_a_y
[1] 1
$ode_a_dydt
[1] 0
$fixed_time_step
[1] 0
$schedule_nsteps
[1] 20
$schedule_eps
[1] 0.02
$schedule_verbose
[1] FALSE
$save_RK45_cache
[1] FALSE
$GSS_tol_abs
[1] 0.001
$vulnerability_curve_ncontrol
[1] 100
$ci_abs_tol
[1] 0.001
$ci_niter
[1] 1000
attr(,"class")
[1] "Control"
$collect_all_auxiliary
[1] FALSE
$birth_rate_x
numeric(0)
$birth_rate_y
[1] 1
$is_variable_birth_rate
[1] FALSE
attr(,"class")
[1] "TF24_Strategy"
$node_schedule_times_default
[1] 0.000000e+00 1.000000e-05 2.000000e-05 3.000000e-05 4.000000e-05
[6] 5.000000e-05 6.000000e-05 7.000000e-05 8.000000e-05 9.525879e-05
[11] 1.105176e-04 1.257764e-04 1.410352e-04 1.562939e-04 1.868115e-04
[16] 2.173291e-04 2.478467e-04 2.783643e-04 3.088818e-04 3.699170e-04
[21] 4.309521e-04 4.919873e-04 5.530225e-04 6.140576e-04 7.361279e-04
[26] 8.581982e-04 9.802686e-04 1.102339e-03 1.224409e-03 1.468550e-03
[31] 1.712690e-03 1.956831e-03 2.200972e-03 2.445112e-03 2.933394e-03
[36] 3.421675e-03 3.909956e-03 4.398237e-03 4.886519e-03 5.863081e-03
[41] 6.839644e-03 7.816206e-03 8.792769e-03 9.769331e-03 1.172246e-02
[46] 1.367558e-02 1.562871e-02 1.758183e-02 1.953496e-02 2.344121e-02
[51] 2.734746e-02 3.125371e-02 3.515996e-02 3.906621e-02 4.687871e-02
[56] 5.469121e-02 6.250371e-02 7.031621e-02 7.812871e-02 9.375371e-02
[61] 1.093787e-01 1.250037e-01 1.406287e-01 1.562537e-01 1.875037e-01
[66] 2.187537e-01 2.500037e-01 2.812537e-01 3.125037e-01 3.750037e-01
[71] 4.375037e-01 5.000037e-01 5.625037e-01 6.250037e-01 7.500037e-01
[76] 8.750037e-01 1.000004e+00 1.125004e+00 1.250004e+00 1.500004e+00
[81] 1.750004e+00 2.000004e+00 2.250004e+00 2.500004e+00 3.000004e+00
[86] 3.500004e+00 4.000004e+00 4.500004e+00 5.000004e+00 6.000004e+00
[91] 7.000004e+00 8.000004e+00 9.000004e+00 1.000000e+01 1.200000e+01
[96] 1.400000e+01 1.600000e+01 1.800000e+01 2.000000e+01 2.200000e+01
[101] 2.400000e+01
$node_schedule_times
$node_schedule_times[[1]]
[1] 0.000000e+00 1.000000e-05 2.000000e-05 3.000000e-05 4.000000e-05
[6] 5.000000e-05 6.000000e-05 7.000000e-05 8.000000e-05 9.525879e-05
[11] 1.105176e-04 1.257764e-04 1.410352e-04 1.562939e-04 1.868115e-04
[16] 2.173291e-04 2.478467e-04 2.783643e-04 3.088818e-04 3.699170e-04
[21] 4.309521e-04 4.919873e-04 5.530225e-04 6.140576e-04 7.361279e-04
[26] 8.581982e-04 9.802686e-04 1.102339e-03 1.224409e-03 1.468550e-03
[31] 1.712690e-03 1.956831e-03 2.200972e-03 2.445112e-03 2.933394e-03
[36] 3.421675e-03 3.909956e-03 4.398237e-03 4.886519e-03 5.863081e-03
[41] 6.839644e-03 7.816206e-03 8.792769e-03 9.769331e-03 1.172246e-02
[46] 1.367558e-02 1.562871e-02 1.758183e-02 1.953496e-02 2.344121e-02
[51] 2.734746e-02 3.125371e-02 3.515996e-02 3.906621e-02 4.687871e-02
[56] 5.469121e-02 6.250371e-02 7.031621e-02 7.812871e-02 9.375371e-02
[61] 1.093787e-01 1.250037e-01 1.406287e-01 1.562537e-01 1.875037e-01
[66] 2.187537e-01 2.500037e-01 2.812537e-01 3.125037e-01 3.750037e-01
[71] 4.375037e-01 5.000037e-01 5.625037e-01 6.250037e-01 7.500037e-01
[76] 8.750037e-01 1.000004e+00 1.125004e+00 1.250004e+00 1.500004e+00
[81] 1.750004e+00 2.000004e+00 2.250004e+00 2.500004e+00 3.000004e+00
[86] 3.500004e+00 4.000004e+00 4.500004e+00 5.000004e+00 6.000004e+00
[91] 7.000004e+00 8.000004e+00 9.000004e+00 1.000000e+01 1.200000e+01
[96] 1.400000e+01 1.600000e+01 1.800000e+01 2.000000e+01 2.200000e+01
[101] 2.400000e+01
$node_schedule_times[[2]]
[1] 0.000000e+00 1.000000e-05 2.000000e-05 3.000000e-05 4.000000e-05
[6] 5.000000e-05 6.000000e-05 7.000000e-05 8.000000e-05 9.525879e-05
[11] 1.105176e-04 1.257764e-04 1.410352e-04 1.562939e-04 1.868115e-04
[16] 2.173291e-04 2.478467e-04 2.783643e-04 3.088818e-04 3.699170e-04
[21] 4.309521e-04 4.919873e-04 5.530225e-04 6.140576e-04 7.361279e-04
[26] 8.581982e-04 9.802686e-04 1.102339e-03 1.224409e-03 1.468550e-03
[31] 1.712690e-03 1.956831e-03 2.200972e-03 2.445112e-03 2.933394e-03
[36] 3.421675e-03 3.909956e-03 4.398237e-03 4.886519e-03 5.863081e-03
[41] 6.839644e-03 7.816206e-03 8.792769e-03 9.769331e-03 1.172246e-02
[46] 1.367558e-02 1.562871e-02 1.758183e-02 1.953496e-02 2.344121e-02
[51] 2.734746e-02 3.125371e-02 3.515996e-02 3.906621e-02 4.687871e-02
[56] 5.469121e-02 6.250371e-02 7.031621e-02 7.812871e-02 9.375371e-02
[61] 1.093787e-01 1.250037e-01 1.406287e-01 1.562537e-01 1.875037e-01
[66] 2.187537e-01 2.500037e-01 2.812537e-01 3.125037e-01 3.750037e-01
[71] 4.375037e-01 5.000037e-01 5.625037e-01 6.250037e-01 7.500037e-01
[76] 8.750037e-01 1.000004e+00 1.125004e+00 1.250004e+00 1.500004e+00
[81] 1.750004e+00 2.000004e+00 2.250004e+00 2.500004e+00 3.000004e+00
[86] 3.500004e+00 4.000004e+00 4.500004e+00 5.000004e+00 6.000004e+00
[91] 7.000004e+00 8.000004e+00 9.000004e+00 1.000000e+01 1.200000e+01
[96] 1.400000e+01 1.600000e+01 1.800000e+01 2.000000e+01 2.200000e+01
[101] 2.400000e+01
$ode_times
numeric(0)
$initial_state
numeric(0)
$n_initial_cohorts
numeric(0)
$initial_node_times
numeric(0)
$initial_patch_density
numeric(0)
$initial_pr_patch_survival
numeric(0)
$initial_time
[1] 0
attr(,"class")
[1] "Parameters<TF24,TF24_Env>" "Parameters" Each line in the size-distribution plot follows one node through time. New nodes enter at times chosen by the node schedule, so a recorded step is not the same as one year. Line transparency represents node density: darker trajectories contribute more individuals to the patch.
results_tidy$species %>%
drop_na() %>%
plot_size_distribution()
Calculate patch totals for each strategy
Integration turns a size distribution into one patch-level value per strategy and time. For example, integrating node density over height gives the total number of individuals per square metre:
data_species_tot <-
results_tidy$species %>% integrate_over_size_distribution()
data_species_tot# A tibble: 202 × 25
step time patch_density species density min_height log_density height
<int> <dbl> <dbl> <chr> <dbl> <dbl> <dbl> <dbl>
1 2 0.00001 0.140 1 0.0000999 0.0351 0.000310 3.51e-6
2 2 0.00001 0.140 2 0.0000999 0.0239 0.000435 2.39e-6
3 3 0.00002 0.140 1 0.000200 0.0351 0.000620 7.01e-6
4 3 0.00002 0.140 2 0.000200 0.0239 0.000870 4.77e-6
5 4 0.00003 0.140 1 0.000300 0.0351 0.000931 1.05e-5
6 4 0.00003 0.140 2 0.000300 0.0239 0.00131 7.16e-6
7 5 0.00004 0.140 1 0.000400 0.0351 0.00124 1.40e-5
8 5 0.00004 0.140 2 0.000400 0.0239 0.00174 9.55e-6
9 6 0.00005 0.140 1 0.000500 0.0351 0.00155 1.75e-5
10 6 0.00005 0.140 2 0.000500 0.0239 0.00218 1.19e-5
# ℹ 192 more rows
# ℹ 17 more variables: mortality <dbl>, fecundity <dbl>, area_heartwood <dbl>,
# mass_heartwood <dbl>, storage <dbl>,
# offspring_produced_survival_weighted <dbl>, competition_effect <dbl>,
# height_inverse <dbl>, net_mass_production_dt <dbl>, root_mass <dbl>,
# opt_psi_stem <dbl>, opt_root_psi <dbl>, transpiration <dbl>, E_up_ <dbl>,
# profit <dbl>, stom_cond_CO2 <dbl>, assimilation <dbl>Plotting that total shows how stem density changes as the patch develops:
data_species_tot %>%
ggplot(aes(time, density, colour=species)) +
geom_line()
expand_state() derives additional quantities, including leaf area and tissue mass, from the states returned by the solver:
results_tidy_expand <- results_tidy %>% expand_state()
results_tidy_expand$species# A tibble: 10,506 × 37
species time step patch_density node density log_density height mortality
<chr> <dbl> <int> <dbl> <int> <dbl> <dbl> <dbl> <dbl>
1 1 0 1 0.140 1 22.3 3.10 0.0351 0.000699
2 1 1e-5 2 0.140 1 22.3 3.10 0.0351 0.000699
3 1 1e-5 2 0.140 2 22.3 3.10 0.0351 0.000699
4 1 2e-5 3 0.140 1 22.3 3.10 0.0351 0.000699
5 1 2e-5 3 0.140 2 22.3 3.10 0.0351 0.000699
6 1 2e-5 3 0.140 3 22.3 3.10 0.0351 0.000699
7 1 3e-5 4 0.140 1 22.3 3.10 0.0351 0.000699
8 1 3e-5 4 0.140 2 22.3 3.10 0.0351 0.000699
9 1 3e-5 4 0.140 3 22.3 3.10 0.0351 0.000699
10 1 3e-5 4 0.140 4 22.3 3.10 0.0351 0.000699
# ℹ 10,496 more rows
# ℹ 28 more variables: fecundity <dbl>, area_heartwood <dbl>,
# mass_heartwood <dbl>, storage <dbl>,
# offspring_produced_survival_weighted <dbl>, competition_effect <dbl>,
# height_inverse <dbl>, net_mass_production_dt <dbl>, root_mass <dbl>,
# opt_psi_stem <dbl>, opt_root_psi <dbl>, transpiration <dbl>, E_up_ <dbl>,
# profit <dbl>, stom_cond_CO2 <dbl>, assimilation <dbl>, area_leaf <dbl>, …To obtain total leaf area, the helper weights individual leaf area by node density and integrates over height at each time:
data_species_tot <-
results_tidy_expand$species %>%
integrate_over_size_distribution()
data_species_tot %>%
ggplot(aes(time, area_leaf, colour = species)) +
geom_line()
In this run, the thin-leaved strategy initially deploys more total leaf area but is later overtaken by the thicker-leaved strategy.
The expanded totals also show how above-ground biomass is divided among leaf, bark, sapwood, and heartwood through time:
v <- c("mass_heartwood", "mass_sapwood", "mass_bark", "mass_leaf")
data_long <-
data_species_tot %>%
select(time,species, one_of(v)) %>%
pivot_longer(cols=starts_with("mass"), names_to = "tissue")
data_long$tissue <- factor(data_long$tissue, levels = v)
species_names <- tibble(species = unique(data_long$species),
species_name = c("Species_1","Species_2"))
data_long %>%
left_join(species_names) %>%
ggplot(aes(time, value, fill=tissue)) +
geom_area() +
labs(x = "Patch age (yr)", y = "Above ground mass (kg/m2)") +
theme_classic() +
xlim(c(0,100)) +
facet_wrap(~species_name)
Interpolate to specific times
The node-introduction schedule sets the times at which new characteristics enter the simulation. It can optionally be refined by the node-spacing algorithm. The recorded steps therefore may not fall at the exact times we want to query. interpolate_to_times() estimates the node states at 1, 5, and 10 years from the existing run:
times <- c(1, 5, 10)
tidy_species_data <- results_tidy$species
tidy_species_new <- interpolate_to_times(tidy_species_data, times)The returned rows contain the state variables at those requested times:
tidy_species_new %>% drop_na()# A tibble: 472 × 24
species node time patch_density density log_density height mortality
<chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1 1 1 0.138 3.57 1.27 2.78 0.0112
2 1 1 5 0.0953 0.404 -0.907 9.49 0.0546
3 1 1 10 0.0298 0.211 -1.55 15.3 0.797
4 1 2 1 0.138 3.57 1.27 2.78 0.0112
5 1 2 5 0.0953 0.404 -0.907 9.48 0.0546
6 1 2 10 0.0298 0.212 -1.55 15.3 0.798
7 1 3 1 0.138 3.57 1.27 2.78 0.0112
8 1 3 5 0.0953 0.404 -0.906 9.48 0.0546
9 1 3 10 0.0298 0.210 -1.56 15.3 0.800
10 1 4 1 0.138 3.57 1.27 2.78 0.0112
# ℹ 462 more rows
# ℹ 16 more variables: fecundity <dbl>, area_heartwood <dbl>,
# mass_heartwood <dbl>, storage <dbl>,
# offspring_produced_survival_weighted <dbl>, competition_effect <dbl>,
# height_inverse <dbl>, net_mass_production_dt <dbl>, root_mass <dbl>,
# opt_psi_stem <dbl>, opt_root_psi <dbl>, transpiration <dbl>, E_up_ <dbl>,
# profit <dbl>, stom_cond_CO2 <dbl>, assimilation <dbl>To check the interpolation, combine the new rows with the recorded trajectories. The red points lie on the original curves, showing agreement at the requested times:
data_combined <-
tidy_species_data %>%
bind_rows(tidy_species_new) %>%
arrange(species, node, time) %>%
filter(node %in% seq(1, 101, by=20))
data_combined_long <-
data_combined %>%
select(node, time, step, density, height, species) %>%
pivot_longer(cols = c("density", "height"))
data_combined_long %>%
ggplot(aes(time, value, group=node,colour=node)) +
geom_line() +
geom_point(data = data_combined_long %>% filter(is.na(step)), col=2) +
scale_y_log10() +
xlim(c(0, 20)) +
facet_grid(name~species, scales="free") +
theme_classic()
Interpolate to specific heights
The complementary question is: what was the state at a particular height? interpolate_to_heights() estimates values across nodes at each recorded time:
heights <- c(1, 5, 10)
tidy_species_data <- results_tidy$species
tidy_species_new <- interpolate_to_heights(tidy_species_data, heights)These rows contain estimates at heights of 1, 5, and 10 metres:
tidy_species_new %>% drop_na()# A tibble: 209 × 24
species time step patch_density density log_density height mortality
<chr> <dbl> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1 0 1 0.140 22.3 3.10 0.0351 0.000699
2 1 0.00001 2 0.140 22.3 3.10 0.0351 0.000699
3 1 0.00002 3 0.140 22.3 3.10 0.0351 0.000699
4 1 0.00003 4 0.140 22.3 3.10 0.0351 0.000699
5 1 0.00004 5 0.140 22.3 3.10 0.0351 0.000699
6 1 0.00005 6 0.140 22.3 3.10 0.0351 0.000700
7 1 0.00006 7 0.140 22.3 3.10 0.0351 0.000700
8 1 0.00007 8 0.140 22.3 3.10 0.0351 0.000700
9 1 0.00008 9 0.140 22.3 3.10 0.0351 0.000700
10 1 0.0000953 10 0.140 22.3 3.10 0.0351 0.000700
# ℹ 199 more rows
# ℹ 16 more variables: fecundity <dbl>, area_heartwood <dbl>,
# mass_heartwood <dbl>, storage <dbl>,
# offspring_produced_survival_weighted <dbl>, competition_effect <dbl>,
# height_inverse <dbl>, net_mass_production_dt <dbl>, root_mass <dbl>,
# opt_psi_stem <dbl>, opt_root_psi <dbl>, transpiration <dbl>, E_up_ <dbl>,
# profit <dbl>, stom_cond_CO2 <dbl>, assimilation <dbl>Again, plotting the estimates as red points against the recorded profiles makes the interpolation easy to assess:
data_combined <-
tidy_species_data %>%
bind_rows(tidy_species_new) %>%
arrange(species, time, height) %>%
filter(step %in% unique(round(seq(min(step, na.rm = TRUE),
max(step, na.rm = TRUE), length.out = 2))))
data_combined_long <-
data_combined %>%
select(node, time, step, height, density, mortality, area_heartwood, species) %>%
pivot_longer(cols = c("density", "mortality", "area_heartwood"))
data_combined_long %>%
filter(!is.na(node)) %>%
ggplot() +
geom_line(aes(height, value, group=step, colour=step)) +
geom_point(data = data_combined_long %>% filter(is.na(node)), aes(x=height, y=value), col="red")+
facet_grid(name~species, scales="free") +
theme_classic()