Two model approaches
The LM3-PPA is a cohort-based vegetation model which simulates vegetation dynamics and biogeochemical processes (Weng et al., 2015). The model is able to link photosynthesis standard models (Farquhar et al., 1980) with tree allometry. In our formulation we retain the original model structure with the standard photosynthesis formulation (i.e. “gs_leuning”) as well as an alternative “p-model” approach. Both model structures operate at different time scales, where the original input has an hourly time step our alternative p-model approach uses a daily time step. Hence, we have two different datasets as driver data (with the lm3ppa p-model input being an aggregate of the high resolution hourly data).
Running the LM3-PPA model with standard photosynthesis
With all data prepared we can run the model using runread_lm3ppa_f(). This function takes the nested data structure and runs the model site by site, returning nested model output results matching the input drivers. In our case only one site will be evaluated.
# print parameter settings
print(lm3ppa_gs_leuning_drivers$params_siml)
#> [[1]]
#> # A tibble: 1 × 12
#> spinup spinupyears recycle firstyeartrend nyeartrend outputhourly outputdaily
#> <lgl> <dbl> <dbl> <dbl> <dbl> <lgl> <lgl>
#> 1 TRUE 250 1 2009 1 TRUE TRUE
#> # … with 5 more variables: do_U_shaped_mortality <lgl>,
#> # update_annualLAImax <lgl>, do_closedN_run <lgl>, method_photosynth <chr>,
#> # method_mortality <chr>
print(head(lm3ppa_gs_leuning_drivers$forcing))
#> [[1]]
#> # A tibble: 8,760 × 13
#> YEAR DOY HOUR PAR Swdown TEMP SoilT RH RAIN WIND PRESSURE
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 2009 1 0 8.17 0.156 0.728 3.16 91.6 0.0000184 3.56 93216.
#> 2 2009 1 1 8.16 0.158 0.780 3.20 91.5 0.0000116 3.37 93189.
#> 3 2009 1 2 8.17 0.162 0.519 3.23 93.2 0.0000116 3.01 93184.
#> 4 2009 1 3 8.15 0.152 0.476 3.28 92.5 0.0000116 3.31 93166.
#> 5 2009 1 4 8.21 0.158 0.336 3.30 92.9 0.0000140 3.23 93143.
#> 6 2009 1 5 8.20 0.161 0.278 3.30 93.8 0.0000140 2.94 93124.
#> 7 2009 1 6 8.18 0.161 0.0966 3.28 95.4 0.0000140 2.98 93114.
#> 8 2009 1 7 8.40 0.164 0.172 3.30 95.4 0.0000211 3.46 93111.
#> 9 2009 1 8 42.1 8.77 0.236 3.33 95.2 0.0000211 3.31 93118.
#> 10 2009 1 9 146. 38.7 0.152 3.41 96.1 0.0000211 3.27 93132.
#> # … with 8,750 more rows, and 2 more variables: aCO2_AW <dbl>, SWC <dbl># run the model
lm3ppa_output_leuning <- runread_lm3ppa_f(
lm3ppa_gs_leuning_drivers,
makecheck = TRUE,
parallel = FALSE
)
# split out the annual data
lm3ppa_gs_leuning_output <- lm3ppa_output_leuning$data[[1]]$output_annual_tilePlotting output
We can now visualize the model output.
# we only have one site so we'll unnest
# the main model output
lm3ppa_gs_leuning_output %>%
ggplot() +
geom_line(aes(x = year, y = GPP)) +
theme_classic()+labs(x = "Year", y = "GPP")
lm3ppa_gs_leuning_output %>%
ggplot() +
geom_line(aes(x = year, y = plantC)) +
theme_classic()+labs(x = "Year", y = "plantC")
Running the LM3-PPA model with P-model photosynthesis
Running the fast P-model implementation.
# print parameter settings
print(lm3ppa_p_model_drivers$params_siml)
#> [[1]]
#> # A tibble: 1 × 12
#> spinup spinupyears recycle firstyeartrend nyeartrend outputhourly outputdaily
#> <lgl> <dbl> <dbl> <dbl> <dbl> <lgl> <lgl>
#> 1 TRUE 250 1 2009 1 TRUE TRUE
#> # … with 5 more variables: do_U_shaped_mortality <lgl>,
#> # update_annualLAImax <lgl>, do_closedN_run <lgl>, method_photosynth <chr>,
#> # method_mortality <chr>
print(head(lm3ppa_p_model_drivers$forcing))
#> [[1]]
#> # A tibble: 365 × 13
#> YEAR DOY HOUR PAR Swdown TEMP SoilT RH RAIN WIND PRESSURE
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 2009 1 11.5 107. 28.0 0.384 3.39 93.7 0.0000166 3.00 93092.
#> 2 2009 2 11.5 102. 24.4 -1.64 2.55 92.5 0.0000232 2.97 93248.
#> 3 2009 3 11.5 197. 46.5 -2.51 1.95 85.1 0.00000371 2.84 93684.
#> 4 2009 4 11.5 139. 34.6 -1.82 2.20 83.5 0.0000130 2.67 93435.
#> 5 2009 5 11.5 154. 36.2 -1.34 2.23 87.8 0.0000223 3.21 93175.
#> 6 2009 6 11.5 99.2 25.9 -0.450 2.75 90.9 0.0000219 3.03 93282.
#> 7 2009 7 11.5 141. 34.8 0.266 3.50 89.7 0.0000136 2.64 93511.
#> 8 2009 8 11.5 169. 44.4 0.504 3.52 86.1 0.0000113 2.68 93443.
#> 9 2009 9 11.5 102. 25.3 0.0869 3.49 90.3 0.0000186 2.74 93447.
#> 10 2009 10 11.5 160. 40.5 -0.404 3.44 90.1 0.0000125 2.17 93633.
#> # … with 355 more rows, and 2 more variables: aCO2_AW <dbl>, SWC <dbl># run the model
lm3ppa_p_model_output <- runread_lm3ppa_f(
lm3ppa_p_model_drivers,
makecheck = TRUE,
parallel = FALSE
)
# split out the annual data for visuals
lm3ppa_p_model_output <- lm3ppa_p_model_output$data[[1]]$output_annual_tilePlotting output
We can now visualize the model output.
# we only have one site so we'll unnest
# the main model output
lm3ppa_p_model_output %>%
ggplot() +
geom_line(aes(x = year, y = GPP)) +
theme_classic()+labs(x = "Year", y = "GPP")
lm3ppa_p_model_output %>%
ggplot() +
geom_line(aes(x = year, y = plantC)) +
theme_classic()+labs(x = "Year", y = "plantC")
Calibrating model parameters
To optimize new parameters based upon driver data and a validation dataset we must first specify an optimization strategy and settings, as well as parameter ranges.
# Mortality as DBH
settings <- list(
method = "bayesiantools",
targetvars = c("gpp"),
timescale = list(targets_obs = "y"),
sitenames = "CH-Lae",
metric = cost_rmse_lm3ppa_gsleuning,
dir_results = "./",
name = "ORG",
control = list(
sampler = "DEzs",
settings = list(
burnin = 10,
iterations = 50
)
),
par = list(
phiRL = list(lower=0.5, upper=5, init=3.5),
LAI_light = list(lower=2, upper=5, init=3.5),
tf_base = list(lower=0.5, upper=1.5, init=1),
par_mort = list(lower=0.1, upper=2, init=1))
)
pars <- calib_sofun(
drivers = lm3ppa_gs_leuning_drivers,
obs = lm3ppa_validation_2,
settings = settings
)