The rsofun package and framework includes two main
models. The pmodel and biomee (which in part
relies on pmodel component). Here we give a short example on how to run
the biomee model on the included demo datasets to
familiarize yourself with both the data structure and the outputs.
Demo data
The package includes two demo datasets to run. These files can be directly loaded into your workspace by typing:
These are real data from the Swiss CH-Lae fluxnet site. We can use
these data to run the model, together with observations of GPP we can
also parameterize biomee parameters.
Two model approaches
BiomeE 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 BiomeE p-model input being an aggregate of the high resolution hourly data).
Running the BiomeE model with standard photosynthesis
With all data prepared we can run the model using
runread_biomee_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
biomee_gs_leuning_drivers$params_siml
#> [[1]]
#> # A tibble: 1 × 11
#> spinup spinupyears recycle firstyeartrend nyeartrend steps_per_day
#> <lgl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 TRUE 250 1 2009 1 24
#> # ℹ 5 more variables: do_U_shaped_mortality <lgl>, update_annualLAImax <lgl>,
#> # do_closedN_run <lgl>, method_photosynth <chr>, method_mortality <chr>
# print forcing
head(biomee_gs_leuning_drivers$forcing)
#> [[1]]
#> # A tibble: 8,760 × 9
#> date hod temp rain vpd ppfd patm wind co2
#> <date> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 2009-01-01 0 0.728 0.0000184 54.4 0.000000319 93216. 3.56 388.
#> 2 2009-01-01 1 0.780 0.0000116 54.9 0.000000322 93189. 3.37 388.
#> 3 2009-01-01 2 0.519 0.0000116 42.9 0.000000330 93184. 3.01 388.
#> 4 2009-01-01 3 0.476 0.0000116 47.3 0.000000311 93166. 3.31 388.
#> 5 2009-01-01 4 0.336 0.0000140 44.3 0.000000322 93143. 3.23 388.
#> 6 2009-01-01 5 0.278 0.0000140 38.4 0.000000329 93124. 2.94 388.
#> 7 2009-01-01 6 0.0966 0.0000140 28.5 0.000000329 93114. 2.98 388.
#> 8 2009-01-01 7 0.172 0.0000211 28.4 0.000000335 93111. 3.46 388.
#> 9 2009-01-01 8 0.236 0.0000211 30.1 0.0000179 93118. 3.31 388.
#> 10 2009-01-01 9 0.152 0.0000211 24.1 0.0000789 93132. 3.27 388.
#> # ℹ 8,750 more rows
set.seed(2023)
# run the model
out <- runread_biomee_f(
biomee_gs_leuning_drivers,
makecheck = TRUE,
parallel = FALSE
)
# split out the annual data
biomee_gs_leuning_output_annual_tile <- out$data[[1]]$output_annual_tile
biomee_gs_leuning_output_annual_cohorts <- out$data[[1]]$output_annual_cohortsPlotting output
We can now visualize the model output.
# we only have one site so we'll unnest
# the main model output
biomee_gs_leuning_output_annual_tile |>
ggplot() +
geom_line(aes(x = year, y = GPP)) +
theme_classic()+labs(x = "Year", y = "GPP")
biomee_gs_leuning_output_annual_tile |>
ggplot() +
geom_line(aes(x = year, y = plantC)) +
theme_classic()+labs(x = "Year", y = "plantC")
biomee_gs_leuning_output_annual_cohorts %>% group_by(cohort,year) %>%
summarise(npp=sum(NPP*density/10000)) %>% mutate(cohort=as.factor(cohort)) %>%
ggplot(aes(x=cohort,y=npp,fill=year)) +
geom_bar(stat="identity") +
theme_classic()+labs(x = "Cohort", y = "NPP")
#> `summarise()` has grouped output by 'cohort'. You can override using the
#> `.groups` argument.
Running the BiomeEP model
Running BiomeE with P-model photosynthesis.
# print parameter settings
biomee_p_model_drivers$params_siml
#> [[1]]
#> # A tibble: 1 × 11
#> spinup spinupyears recycle firstyeartrend nyeartrend steps_per_day
#> <lgl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 TRUE 250 1 2009 1 1
#> # ℹ 5 more variables: do_U_shaped_mortality <lgl>, update_annualLAImax <lgl>,
#> # do_closedN_run <lgl>, method_photosynth <chr>, method_mortality <chr>
# print forcing for P-model
head(biomee_p_model_drivers$forcing)
#> [[1]]
#> # A tibble: 365 × 9
#> date hod temp rain vpd ppfd patm wind co2
#> <date> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 2009-01-01 11.5 0.384 0.0000166 39.5 0.0000571 93092. 3.00 388.
#> 2 2009-01-02 11.5 -1.64 0.0000232 40.5 0.0000499 93248. 2.97 388.
#> 3 2009-01-03 11.5 -2.51 0.00000371 75.9 0.0000949 93684. 2.84 388.
#> 4 2009-01-04 11.5 -1.82 0.0000130 88.0 0.0000706 93435. 2.67 388.
#> 5 2009-01-05 11.5 -1.34 0.0000223 67.8 0.0000739 93175. 3.21 388.
#> 6 2009-01-06 11.5 -0.450 0.0000219 54.0 0.0000528 93282. 3.03 388.
#> 7 2009-01-07 11.5 0.266 0.0000136 64.1 0.0000711 93511. 2.64 388.
#> 8 2009-01-08 11.5 0.504 0.0000113 88.2 0.0000905 93443. 2.68 388.
#> 9 2009-01-09 11.5 0.0869 0.0000186 59.9 0.0000516 93447. 2.74 388.
#> 10 2009-01-10 11.5 -0.404 0.0000125 58.6 0.0000827 93633. 2.17 388.
#> # ℹ 355 more rows
# run the model
out <- runread_biomee_f(
biomee_p_model_drivers,
makecheck = TRUE,
parallel = FALSE
)
# split out the annual data for visuals
biomee_p_model_output_annual_tile <- out$data[[1]]$output_annual_tile
biomee_p_model_output_annual_cohorts <- out$data[[1]]$output_annual_cohortsPlotting output
We can now visualize the model output.
# we only have one site so we'll unnest
# the main model output
biomee_p_model_output_annual_tile %>%
ggplot() +
geom_line(aes(x = year, y = GPP)) +
theme_classic() +
labs(x = "Year", y = "GPP")
biomee_p_model_output_annual_tile %>%
ggplot() +
geom_line(aes(x = year, y = plantC)) +
theme_classic() +
labs(x = "Year", y = "plantC")
biomee_p_model_output_annual_cohorts %>% group_by(cohort,year) %>%
summarise(npp=sum(NPP*density/10000)) %>% mutate(cohort=as.factor(cohort)) %>%
ggplot(aes(x=cohort,y=npp,fill=year)) +
geom_bar(stat="identity") +
theme_classic()+labs(x = "Cohort", y = "NPP")
#> `summarise()` has grouped output by 'cohort'. You can override using the
#> `.groups` argument.
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. In this example, we use as cost the root mean
squared error (RMSE) between simulated and observed targets (GPP, LAI,
Density and Biomass) and we minimize it using the GenSA
optimizer.
# Mortality as DBH
settings <- list(
method = "GenSA",
metric = cost_rmse_biomee,
control = list(
maxit = 10
),
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))
)
# Using BiomeEP (with P-model for photosynthesis)
pars <- calib_sofun(
drivers = biomee_p_model_drivers,
obs = biomee_validation,
settings = settings
)Using the calibrated parameter values, we can run again the BiomeE simulations.
# replace parameter values by calibration output
drivers <- biomee_p_model_drivers
drivers$params_species[[1]]$phiRL[] <- pars$par[1]
drivers$params_species[[1]]$LAI_light[] <- pars$par[2]
drivers$params_tile[[1]]$tf_base <- pars$par[3]
drivers$params_tile[[1]]$par_mort <- pars$par[4]
# run the model with new parameter values
biomee_p_model_output_calib <- runread_biomee_f(
drivers,
makecheck = TRUE,
parallel = FALSE
)
# split out the annual data
biomee_p_model_output_calib <- biomee_p_model_output_calib$data[[1]]$output_annual_tileFinally, we can visually compare the two model runs. We plot the original model run in black and the run using the calibrated parameters in grey. The dashed line represents the validation data, i.e. the observed GPP and Biomass.
# unnest model output for our single site
ggplot() +
geom_line(data = biomee_p_model_output$data[[1]]$output_annual_tile,
aes(x = year, y = GPP)) +
geom_line(data = biomee_p_model_output_calib,
aes(x = year, y = GPP),
color = "grey50") +
theme_classic() +
labs(x = "Year", y = "GPP")
ggplot() +
geom_line(data = biomee_p_model_output$data[[1]]$output_annual_tile,
aes(x = year, y = plantC)) +
geom_line(data = biomee_p_model_output_calib,
aes(x = year, y = plantC),
color = "grey50") +
theme_classic() +
labs(x = "Year", y = "plantC")