Title: | The P-Model and BiomeE Modelling Framework |
---|---|
Description: | Implements the Simulating Optimal FUNctioning framework for site-scale simulations of ecosystem processes, including model calibration. It contains 'Fortran 90' modules for the P-model (Stocker et al. (2020) <doi:10.5194/gmd-13-1545-2020>), SPLASH (Davis et al. (2017) <doi:10.5194/gmd-10-689-2017>) and BiomeE (Weng et al. (2015) <doi:10.5194/bg-12-2655-2015>). |
Authors: | Benjamin Stocker [aut, cre] , Koen Hufkens [aut] , Josefa Arán Paredes [aut] , Laura Marqués [ctb] , Mayeul Marcadella [ctb] , Ensheng Weng [ctb] , Fabian Bernhard [aut] , Geocomputation and Earth Observation, University of Bern [cph, fnd] |
Maintainer: | Benjamin Stocker <[email protected]> |
License: | GPL-3 |
Version: | 5.0.0.9000 |
Built: | 2024-12-12 15:29:44 UTC |
Source: | https://github.com/geco-bern/rsofun |
Small dataset representing the driver to run the BiomeE-model at the CH-LAE site
using the Leuning photosynthesis specification (and half-hourly time step)
It can also be used together with leaf trait data from CH-LAE (biomee_validation
)
to optimize model parameters.
biomee_gs_leuning_drivers
biomee_gs_leuning_drivers
A tibble of driver data.
Site name
Simulation parameters as a data.frame, including the following data:
Flag indicating whether this simulation does spin-up.
Number of spin-up years.
Number of first N years of forcing data.frame that are recycled for spin-up.
Year of first transient year (AD) (optional). Is only used to set years in output data frames. Defaults to 0 if not provided.
Number of transient years (optional). Determines the length of simulation output after spin-up. Defaults to number of years contained in the forcing data. (If longer than forcing data, last year of forcing is repeated until the end (spin-down).)
Time resolution of the forcing (day-1).
Flag indicating whether U-shaped mortality is used.
Flag indicating whether updating LAImax according to mineral N in soil.
Flag indicating whether doing N closed runs to recover N balance enforcing 0.2 kg N m-2 in the inorganic N pool.
String specifying the method of photosynthesis used in the model, either "pmodel" or "gs_leuning".document()
String indicating the type of mortality in the model. One of the following: "dbh" is size-dependent mortality, "const_selfthin" is constant self thinning (in development), "cstarvation" is carbon starvation, and "growthrate" is growth rate dependent mortality.
Site meta info in a data.frame. This data structure can be freely used for documenting the dataset, but must include at least the following data:
Longitude of the site location in degrees east.
Latitude of the site location in degrees north.
Elevation of the site location, in meters above sea level.
Forcing data.frame used as input
Photosynthetic photon flux density (mol s-1 m-2)
Air temperature (deg C)
Vapor pressure deficit (Pa)
Precipitation (kgH2O m-2 s-1 == mm s-1)
Wind velocity (m s-1)
Atmospheric pressure (pa)
CO2 atmospheric concentration (ppm)
Tile-level model parameters, into a single row data.frame, including the following data:
Integer indicating the type of soil: Sand = 1, LoamySand = 2, SandyLoam = 3, SiltLoam = 4, FrittedClay = 5, Loam = 6, Clay = 7.
Field capacity (vol/vol). Water remaining in a soil after it has been thoroughly saturated and allowed to drain freely.
Wilting point (vol/vol). Water content of a soil at which plants wilt and fail to recover.
Fast soil C decomposition rate (year).
Slow soil C decomposition rate (year).
Mineral Nitrogen turnover rate (year).
Ratio of C and N returned to litters from microbes.
N loss rate through runoff (organic and mineral) (year).
Minimum LMA, leaf mass per unit area, kg C m.
Fraction of fast turnover carbon in fine biomass.
Fraction of fast turnover carbon in wood biomass.
Growth respiration factor.
Fraction of the carbon retained after leaf drop.
Retranslocation coefficient of nitrogen.
Coefficient for setting up initial sapwood.
Re-fill of N for sapwood.
Calibratable scalar for respiration, used to increase LUE levels.
Canopy mortality parameter.
Parameter for understory mortality.
A data.frame containing species-specific model parameters, with one species per row, including the following data:
Integer set to 0 for grasses and 1 for trees.
Integer set to 0 for deciduous and 1 for evergreen.
Integer indicating the type of plant according to photosynthesis: 0 for C3; 1 for C4
Fine root turnover rate (year).
Material density of fine roots (kg C m).
Radius of the fine roots, in m.
e-folding parameter of root vertical distribution, in m.
Fine root water conductivity (mol m
s
MPa
).
Characteristic leaf size.
Max RuBisCo rate, in mol m s
.
Annual productivity per unit area at full sun (kg C
m year
).
Wet leaf photosynthesis down-regulation.
Factor of stomatal conductance.
Photosynthesis efficiency.
Leaf respiration coefficient, in year.
Leaf respiration coefficient per unit N.
Sapwood respiration rate, in kg C m year
.
Fine root respiration rate, kg C kg C
year
.
Critical temperature triggerng offset of phenology, in Kelvin.
Critical temperature triggerng onset of phenology, in Kelvin.
Critical value of GDD5 for turning ON growth season.
Critical soil moisture for phenology onset.
Critical soil moisture for phenology offset.
Initial size of seedlings, in kg C per individual.
Basal leaf N per unit area, in kg N m.
Maximum crown LAI (leaf area index).
Reference N fixation rate (kg N kg C root).
Carbon cost of N fixation (kg C kg N).
Ratio of sapwood area to leaf area.
Canopy tree mortality rate (year).
Understory tree mortality rate (year).
Age at which trees can reproduce (years).
Fraction of G_SF to G_F.
Multiplier for NSNmax as sum of potential bl and br.
Leaf mass per unit area (kg C m).
Wood density (kg C m).
Coefficient for allometry (biomass = alphaBM * DBH ** thetaBM).
Coefficient for allometry (biomass = alphaBM * DBH ** thetaBM).
Quantum yield efficiency ,
in mol mol
.
Ratio of fine root to leaf area.
Maximum LAI limited by light.
A data.frame of initial cohort specifications, including the following data:
Index of a species described in param_species.
Initial individual density, in individuals per
m.
Initial biomass of leaf, in kg C per individual.
Initial biomass of fine root, in kg C per individual.
Initial biomass of sapwood, in kg C per individual.
Initial biomass of heartwood, in kg C per individual.
Initial biomass of seed, in kg C per individual.
Initial non-structural biomass, in kg C per individual.
A data.frame of initial soil pools, including the following data:
Initial fast soil carbon, in kg C m.
Initial slow soil carbon, in kg C m.
Mineral nitrogen pool, in kg N m.
Annual nitrogen input to soil N pool, in kg N m
year
.
Example output dataset from a BiomeE-model run (gs_leuning)
See run_biomee_f_bysite
for a detailed
description of the outputs.
biomee_gs_leuning_output
biomee_gs_leuning_output
An object of class tbl_df
(inherits from tbl
, data.frame
) with 1 rows and 2 columns.
Small dataset representing the driver to run the BiomeE-model at the CH-LAE site
using the P-model photosynthesis specification (and daily time step).
It can also be used together with leaf trait data from CH-LAE (biomee_validation
)
to optimize model parameters.
biomee_p_model_drivers
biomee_p_model_drivers
Example output dataset from a BiomeE-model run (p-model)
See run_biomee_f_bysite
for a detailed
description of the outputs.
biomee_p_model_output
biomee_p_model_output
An object of class tbl_df
(inherits from tbl
, data.frame
) with 1 rows and 2 columns.
Small example dataset of target observations (leaf trait data) at the CH-LAE site
to optimize model parameters with the function calib_sofun
biomee_validation
biomee_validation
A tibble of validation data:
site name
validation data
Lukas Hörtnagl, Werner Eugster, Nina Buchmann, Eugenie Paul-Limoges, Sophia Etzold, Matthias Haeni, Peter Pluess, Thomas Baur (2004-2014) FLUXNET2015 CH-Lae Laegern, Dataset. https://doi.org/10.18140/FLX/1440134
This is the main function that handles the calibration of SOFUN model parameters.
calib_sofun(drivers, obs, settings, optim_out = TRUE, ...)
calib_sofun(drivers, obs, settings, optim_out = TRUE, ...)
drivers |
A data frame with driver data. See |
obs |
A data frame containing observational data used for model
calibration. See |
settings |
A list containing model calibration settings. See the 'P-model usage' vignette for more information and examples.
|
optim_out |
A logical indicating whether the function returns the raw output of the optimization functions (defaults to TRUE). |
... |
Optional arguments passed on to the cost function specified as
|
A named list containing the calibrated parameter vector 'par' and the output object from the optimization 'mod'. For more details on this output and how to evaluate it, see runMCMC (also this post) and GenSA.
# Fix model parameters that won't be calibrated params_fix <- list( kphio_par_a = 0, kphio_par_b = 1.0, soilm_thetastar = 0.6*240, soilm_betao = 0.01, beta_unitcostratio = 146, rd_to_vcmax = 0.014, tau_acclim = 30, kc_jmax = 0.41 ) # Define calibration settings settings <- list( method = "BayesianTools", par = list( kphio = list(lower=0.04, upper=0.09, init=0.05), err_gpp = list(lower = 0.01, upper = 4, init = 2) ), metric = rsofun::cost_likelihood_pmodel, control = list( sampler = "DEzs", settings = list( nrChains = 1, burnin = 0, iterations = 50 # kept artificially low ) ) ) # Run the calibration for GPP data calib_output <- rsofun::calib_sofun( drivers = rsofun::p_model_drivers, obs = rsofun::p_model_validation, settings = settings, # extra arguments for the cost function par_fixed = params_fix, targets = c("gpp") )
# Fix model parameters that won't be calibrated params_fix <- list( kphio_par_a = 0, kphio_par_b = 1.0, soilm_thetastar = 0.6*240, soilm_betao = 0.01, beta_unitcostratio = 146, rd_to_vcmax = 0.014, tau_acclim = 30, kc_jmax = 0.41 ) # Define calibration settings settings <- list( method = "BayesianTools", par = list( kphio = list(lower=0.04, upper=0.09, init=0.05), err_gpp = list(lower = 0.01, upper = 4, init = 2) ), metric = rsofun::cost_likelihood_pmodel, control = list( sampler = "DEzs", settings = list( nrChains = 1, burnin = 0, iterations = 50 # kept artificially low ) ) ) # Run the calibration for GPP data calib_output <- rsofun::calib_sofun( drivers = rsofun::p_model_drivers, obs = rsofun::p_model_validation, settings = settings, # extra arguments for the cost function par_fixed = params_fix, targets = c("gpp") )
Cost function for parameter calibration, which computes the log-likelihood for the biomee model fitting several target variables for a given set of parameters.
cost_likelihood_biomee(par, obs, drivers, targets)
cost_likelihood_biomee(par, obs, drivers, targets)
par |
A vector containing parameter values for |
obs |
A nested data frame of observations, following the structure of |
drivers |
A nested data frame of driver data, for example |
targets |
A character vector indicating the target variables for which the
optimization will be done. This should be a subset of |
The cost function performs a BiomeE model run for the value of
par
given as argument. The likelihood is calculated assuming that the
predicted targets are independent, normally distributed and centered on the observations.
The optimization
should be run using BayesianTools
, so the likelihood is maximized.
The log-likelihood of the simulated targets by the biomee model versus the observed targets.
# Compute the likelihood for a set of # BiomeE model parameter values # and the example data cost_likelihood_biomee( par = c(3.5, 3.5, 1, 1, # model params 0.5), # err_GPP obs = biomee_validation, drivers = biomee_gs_leuning_drivers, targets = c("GPP") )
# Compute the likelihood for a set of # BiomeE model parameter values # and the example data cost_likelihood_biomee( par = c(3.5, 3.5, 1, 1, # model params 0.5), # err_GPP obs = biomee_validation, drivers = biomee_gs_leuning_drivers, targets = c("GPP") )
The cost function performs a P-model run for the input drivers and model parameter values, and computes the outcome's normal log-likelihood centered at the input observed values and with standard deviation given as an input parameter (calibratable).
cost_likelihood_pmodel( par, obs, drivers, targets, par_fixed = NULL, parallel = FALSE, ncores = 2 )
cost_likelihood_pmodel( par, obs, drivers, targets, par_fixed = NULL, parallel = FALSE, ncores = 2 )
par |
A vector of values for the parameters to be calibrated, including
a subset of model parameters (described in |
obs |
A nested data.frame of observations, with columns |
drivers |
A nested data.frame of driver data. See |
targets |
A character vector indicating the target variables for which the
optimization will be done and the RMSE computed. This string must be a column
name of the |
par_fixed |
A named list of model parameter values to keep fixed during the
calibration. These should complement the input |
parallel |
A logical specifying whether simulations are to be parallelised
(sending data from a certain number of sites to each core). Defaults to
|
ncores |
An integer specifying the number of cores used for parallel computing. Defaults to 2. |
To run the P-model, all model parameters must be given. The cost
function uses arguments par
and par_fixed
such that, in the
calibration routine, par
can be updated by the optimizer and
par_fixed
are kept unchanged throughout calibration.
If the validation data contains a "date" column (fluxes), the simulated target time series is compared to the observed values on those same dates (e.g. for GPP). Otherwise, there should only be one observed value per site (leaf traits), and the outputs (averaged over the growing season, weighted by predicted GPP) will be compared to this single value representative of the site (e.g. Vcmax25). As an exception, when the date of a trait measurement is available, it will be compared to the trait value predicted on that date.
The log-likelihood of the observed target values, assuming that they are independent, normally distributed and centered on the predictions made by the P-model run with standard deviation given as input (via 'par' because the error terms are estimated through the calibration with 'BayesianTools', as shown in the "Parameter calibration and cost functions" vignette).
# Compute the likelihood for a set of # model parameter values involved in the # temperature dependence of kphio # and example data cost_likelihood_pmodel( par = c(0.05, -0.01, 1, # model parameters 2), # err_gpp obs = p_model_validation, drivers = p_model_drivers, targets = c('gpp'), par_fixed = list( soilm_thetastar = 0.6 * 240, # old setup with soil moisture stress soilm_betao = 0.0, beta_unitcostratio = 146.0, rd_to_vcmax = 0.014, # from Atkin et al. 2015 for C3 herbaceous tau_acclim = 30.0, kc_jmax = 0.41 ) )
# Compute the likelihood for a set of # model parameter values involved in the # temperature dependence of kphio # and example data cost_likelihood_pmodel( par = c(0.05, -0.01, 1, # model parameters 2), # err_gpp obs = p_model_validation, drivers = p_model_drivers, targets = c('gpp'), par_fixed = list( soilm_thetastar = 0.6 * 240, # old setup with soil moisture stress soilm_betao = 0.0, beta_unitcostratio = 146.0, rd_to_vcmax = 0.014, # from Atkin et al. 2015 for C3 herbaceous tau_acclim = 30.0, kc_jmax = 0.41 ) )
Cost function for parameter calibration, which computes the root mean squared
error (RMSE) between BiomeE simulations (using the input set of parameters)
and observed target variables.
Cost function for parameter calibration, which
computes the RMSE for the biomee model fitting target variables
'GPP','LAI','Density'
and 'Biomass'
for a given set of parameters.
cost_rmse_biomee(par, obs, drivers)
cost_rmse_biomee(par, obs, drivers)
par |
A vector containing parameter values for |
obs |
A nested data frame of observations, following the structure of |
drivers |
A nested data frame of driver data, for example |
The root mean squared error (RMSE) between the observed and simulated
values of 'GPP','LAI','Density'
and 'Biomass'
(all variables
have the same weight). Relative errors (difference divided by observed values) are used
instead of absolute errors.
The cost function performs a BiomeE model run for parameter values
par
and model drivers drivers
given as arguments, producing the
simulated values used to compute the RMSE.
# Compute RMSE for a set of # model parameter values # and example data cost_rmse_biomee( par = c(3.5, 3.5, 1, 1), obs = biomee_validation, drivers = biomee_gs_leuning_drivers )
# Compute RMSE for a set of # model parameter values # and example data cost_rmse_biomee( par = c(3.5, 3.5, 1, 1), obs = biomee_validation, drivers = biomee_gs_leuning_drivers )
The cost function performs a P-model run for the input drivers and parameter values, and compares the output to observations of various targets by computing the root mean squared error (RMSE).
cost_rmse_pmodel( par, obs, drivers, targets, par_fixed = NULL, target_weights = NULL, parallel = FALSE, ncores = 2 )
cost_rmse_pmodel( par, obs, drivers, targets, par_fixed = NULL, target_weights = NULL, parallel = FALSE, ncores = 2 )
par |
A vector of values for the parameters to be calibrated (a subset of
those described in |
obs |
A nested data.frame of observations, with columns |
drivers |
A nested data.frame of driver data. See |
targets |
A character vector indicating the target variables for which the
optimization will be done and the RMSE computed. This string must be a column
name of the |
par_fixed |
A named list of model parameter values to keep fixed during the
calibration. These should complement the input |
target_weights |
A vector of weights to be used in the computation of
the RMSE if using several targets. By default ( |
parallel |
A logical specifying whether simulations are to be parallelised
(sending data from a certain number of sites to each core). Defaults to
|
ncores |
An integer specifying the number of cores used for parallel computing. Defaults to 2. |
To run the P-model, all model parameters must be given. The cost
function uses arguments par
and par_fixed
such that, in the
calibration routine, par
can be updated by the optimizer and
par_fixed
are kept unchanged throughout calibration.
If the validation data contains a "date" column (fluxes), the simulated target time series is compared to the observed values on those same dates (e.g. for GPP). Otherwise, there should only be one observed value per site (leaf traits), and the outputs (averaged over the growing season, weighted by predicted GPP) will be compared to this single value representative of the site (e.g. Vcmax25). As an exception, when the date of a trait measurement is available, it will be compared to the trait value predicted on that date.
The root mean squared error (RMSE) between observed values and P-model predictions. The RMSE is computed for each target separately and then aggregated (mean or weighted average).
# Compute RMSE for a set # of model parameter values # and example data cost_rmse_pmodel( par = c(0.05, -0.01, 0.5), # kphio related parameters obs = p_model_validation, drivers = p_model_drivers, targets = c('gpp'), par_fixed = list( soilm_thetastar = 0.6 * 240, # old setup with soil moisture stress soilm_betao = 0.0, beta_unitcostratio = 146.0, rd_to_vcmax = 0.014, # from Atkin et al. 2015 for C3 herbaceous tau_acclim = 30.0, kc_jmax = 0.41 ) )
# Compute RMSE for a set # of model parameter values # and example data cost_rmse_pmodel( par = c(0.05, -0.01, 0.5), # kphio related parameters obs = p_model_validation, drivers = p_model_drivers, targets = c('gpp'), par_fixed = list( soilm_thetastar = 0.6 * 240, # old setup with soil moisture stress soilm_betao = 0.0, beta_unitcostratio = 146.0, rd_to_vcmax = 0.014, # from Atkin et al. 2015 for C3 herbaceous tau_acclim = 30.0, kc_jmax = 0.41 ) )
Creates a tibble with rows for each date from 'yrstart'
to 'yrend'
in 'yyyy-mm-dd'
format. Intervals of dates are specified by argument
'freq'
.
ddf <- init_dates_dataframe(2000, 2003, startmoy=1, startdoy=1,
freq="days", endmoy=12, enddom=31, noleap=FALSE)
init_dates_dataframe( yrstart, yrend, startmoy = 1, startdoy = 1, freq = "days", endmoy = 12, enddom = 31, noleap = FALSE )
init_dates_dataframe( yrstart, yrend, startmoy = 1, startdoy = 1, freq = "days", endmoy = 12, enddom = 31, noleap = FALSE )
yrstart |
An integer defining the start year of dates covered by the dataframe. |
yrend |
An integer defining the end year of dates covered by the dataframe. |
startmoy |
An integer defining the start month-of-year of dates covered by the dataframe. Defaults to 1. |
startdoy |
An integer defining the start day-of-year of dates covered by the dataframe. Defaults to 1. |
freq |
A character string specifying the time steps of dates
(in rows). Defaults to |
endmoy |
An integer defining the end month-of-year of dates covered by the dataframe. Defaults to 12. |
enddom |
An integer defining the end day-of-year of dates covered by the dataframe. Defaults to 31. |
noleap |
Whether leap years are ignored, that is, whether the 29 |
A tibble with dates.
Small dataset representing the driver to run the P-model at the FR-Pue site.
It can also be used together with daily GPP flux time series data from CH-LAE
(p_model_validation
) to optimize model parameters.
To optimize model parameters to leaf traits data use the datasets p_model_drivers_vcmax25
and p_model_validation_vcmax25
.
p_model_drivers
p_model_drivers
A tibble of driver data:
A character string containing the site name.
A tibble of a time series of forcing climate data, including the following data:
Date of the observation in YYYY-MM-DD format.
Daytime average air temperature in C.
Daytime average vapour pressure deficit in Pa.
Photosynthetic photon flux density (PPFD) in
mol m s
. If all values are NA, it indicates that
PPFD should be calculated by the SPLASH model.
Net radiation in W m. This is currently
ignored as a model forcing.
Atmospheric pressure in Pa.
Snow in water equivalents mm s.
Rain as precipitation in liquid form in mm s.
Daily minimum air temperature in C.
Daily maximum air temperature in C.
Fraction of photosynthetic active radiation (fAPAR), taking values between 0 and 1.
Atmospheric CO concentration.
Cloud coverage in %. This is only used when either PPFD or net radiation are not prescribed.
A tibble of simulation parameters, including the following data:
A logical value indicating whether this simulation does spin-up.
Number of spin-up years.
Number of first N years of forcing data.frame that are recycled for spin-up.
An integer indicating the output periodicity.
A logical value, TRUE
if evergreen tree.
A logical value, TRUE
if evergreen tree and N-fixing.
A logical value, TRUE
if deciduous tree.
A logical value, TRUE
if deciduous tree and N-fixing.
A logical value, TRUE
if grass with C3 photosynthetic pathway.
A logical value, TRUE
if grass with C3 photosynthetic
pathway and N-fixing.
A logical value, TRUE
if grass with C4 photosynthetic pathway.
A tibble containing site meta information. This data structure can be freely used for documenting the dataset, but must include at least the following data:
Longitude of the site location in degrees east.
Latitude of the site location in degrees north.
Elevation of the site location, in meters above sea level.
A numeric value for the rooting zone water holding capacity (in mm)
Pastorello, G., Trotta, C., Canfora, E. et al. The FLUXNET2015 dataset and the ONEFlux processing pipeline for eddy covariance data. Sci Data 7, 225 (2020). https://doi.org/10.1038/s41597-020-0534-3
University of East Anglia Climatic Research Unit; Harris, I.C.; Jones, P.D.; Osborn, T. (2021): CRU TS4.05: Climatic Research Unit (CRU) Time-Series (TS) version 4.05 of high-resolution gridded data of month-by-month variation in climate (Jan. 1901- Dec. 2020). NERC EDS Centre for Environmental Data Analysis, date of citation. https://catalogue.ceda.ac.uk/uuid/c26a65020a5e4b80b20018f148556681
Weedon, G. P., G. Balsamo, N. Bellouin,S. Gomes, M. J. Best, and P. Viterbo(2014), The WFDEI meteorologicalforcing data set: WATCH Forcing Datamethodology applied to ERA-Interimreanalysis data, Water Resour. Res.,50,7505–7514, doi:10.1002/2014WR015638.
Fick, S.E. and R.J. Hijmans, 2017. WorldClim 2: new 1km spatial resolution climate surfaces for global land areas. International Journal of Climatology 37 (12): 4302-4315.
Small dataset representing the driver to run the P-model at four separate sites.
It can also be used together with leaf traits data from these four sites
(p_model_validation_vcmax25
) to optimize model parameters.
To optimize model parameters to GPP flux data use the datasets p_model_drivers
and p_model_validation
.
p_model_drivers_vcmax25
p_model_drivers_vcmax25
See p_model_drivers
Atkin, O. K., Bloomfield, K. J., Reich, P. B., Tjoelker, M. G., Asner, G. P., Bonal, D., et al. (2015). Global variability in leaf respiration in relation to climate, plant functional types and leaf traits. New Phytol. 206 (2), 614–636. doi:10.1111/nph.13253
University of East Anglia Climatic Research Unit; Harris, I.C.; Jones, P.D.; Osborn, T. (2021): CRU TS4.05: Climatic Research Unit (CRU) Time-Series (TS) version 4.05 of high-resolution gridded data of month-by-month variation in climate (Jan. 1901- Dec. 2020). NERC EDS Centre for Environmental Data Analysis, date of citation. https://catalogue.ceda.ac.uk/uuid/c26a65020a5e4b80b20018f148556681
Weedon, G. P., G. Balsamo, N. Bellouin,S. Gomes, M. J. Best, and P. Viterbo(2014), The WFDEI meteorologicalforcing data set: WATCH Forcing Datamethodology applied to ERA-Interimreanalysis data, Water Resour. Res.,50,7505–7514, doi:10.1002/2014WR015638.
Fick, S.E. and R.J. Hijmans, 2017. WorldClim 2: new 1km spatial resolution climate surfaces for global land areas. International Journal of Climatology 37 (12): 4302-4315.
C.D. Keeling, R.B. Bacastow, A.E. Bainbridge, C.A. Ekdahl, P.R. Guenther, and L.S. Waterman, (1976), Atmospheric carbon dioxide variations at Mauna Loa Observatory, Hawaii, Tellus, vol. 28, 538-551
Example output dataset from a p-model run using p_model_drivers
See run_pmodel_f_bysite
for a detailed
description of the outputs.
p_model_output
p_model_output
An object of class tbl_df
(inherits from tbl
, data.frame
) with 1 rows and 3 columns.
Example output dataset from a p-model run using p_model_drivers_vcmax25
See run_pmodel_f_bysite
for a detailed
description of the outputs.
p_model_output_vcmax25
p_model_output_vcmax25
An object of class tbl_df
(inherits from tbl
, data.frame
) with 4 rows and 3 columns.
Small example dataset of target observations (daily GPP flux data) to optimize
model parameters with the function calib_sofun
p_model_validation
p_model_validation
A tibble of validation data:
A character string containing the site name (e.g. 'FR-Pue').
A tibble [ 2,920 x 3 ] with time series for the following variables:
Date vector with format YYYY-MM-DD.
The observed Gross Primary Productivity (GPP) for each time stamp
(in gC m d
).
The uncertainty of the GPP (in gC m d
).
Pastorello, G., Trotta, C., Canfora, E. et al. The FLUXNET2015 dataset and the ONEFlux processing pipeline for eddy covariance data. Sci Data 7, 225 (2020). https://doi.org/10.1038/s41597-020-0534-3
require(ggplot2); require(tidyr) p_model_validation %>% tidyr::unnest(data)
require(ggplot2); require(tidyr) p_model_validation %>% tidyr::unnest(data)
Small example dataset of target observations (leaf trait data) to optimize
model parameters with the function calib_sofun
p_model_validation_vcmax25
p_model_validation_vcmax25
A tibble of validation data:
A character string containing the site names (e.g. 'Reichetal_Colorado').
A tibble [ 1 x 2 ] with observations for the following variables:
The observed maximum rate of carboxylation (Vcmax), normalised
to 25 C (in mol C m
d
), aggregated over different plant species
in each site.
The uncertainty of the Vcmax25 (in mol C m d
),
calculated as the standard deviation among Vcmax25 observations for
several species per site or as the total standard deviation across sites for
single-plant-species sites.
Atkin, O. K., Bloomfield, K. J., Reich, P. B., Tjoelker, M. G., Asner, G. P., Bonal, D., et al. (2015). Global variability in leaf respiration in relation to climate, plant functional types and leaf traits. New Phytol. 206 (2), 614–636. doi:10.1111/nph.13253
require(ggplot2); require(tidyr) p_model_validation_vcmax25 %>% tidyr::unnest(data)
require(ggplot2); require(tidyr) p_model_validation_vcmax25 %>% tidyr::unnest(data)
Run BiomeE Fortran model on single site.
run_biomee_f_bysite( sitename, params_siml, site_info, forcing, params_tile, params_species, init_cohort, init_soil, makecheck = TRUE )
run_biomee_f_bysite( sitename, params_siml, site_info, forcing, params_tile, params_species, init_cohort, init_soil, makecheck = TRUE )
sitename |
Site name. |
params_siml |
Simulation parameters. |
site_info |
Site meta info in a data.frame. |
forcing |
A data frame of forcing climate data, used as input. |
params_tile |
Tile-level model parameters, into a single row data.frame. |
params_species |
A data.frame containing species-specific model parameters,
with one species per row. See examples |
init_cohort |
A data.frame of initial cohort specifications. |
init_soil |
A data.frame of initial soil pools. |
makecheck |
A logical specifying whether checks are performed
to verify forcings and model parameters. For further specifications of above inputs and examples see |
Model output is provided as a list, with elements:
output_hourly_tile
A data.frame containing hourly predictions .
Year of the simulation.
Day of the year.
Hour of the day.
Radiation, in W m.
Air temperature, in Kelvin.
Precipitation, in mm m.
Gross primary production (kg C m hour
).
Plant respiration (kg C m hour
).
Transpiration (mm m).
Evaporation (mm m).
Water runoff (mm m).
Soil water content in root zone (kg m).
Volumetric soil water content for each layer (vol/vol).
Field capacity (vol/vol).
Wilting point (vol/vol).
output_daily_tile
A data.frame with daily outputs at a tile level.
Year of the simulation.
Day of the year.
Air temperature (Kelvin).
Precipitation (mm m).
Soil water content in root zone (kg m).
Transpiration (mm m).
Evaporation (mm m).
Water runoff (mm m).
Volumetric soil water content for layer 1.
Volumetric soil water content for layer 2.
Volumetric soil water content for layer 3.
Leaf area index (m/m
).
Gross primary production (kg C m day
).
Plant autotrophic respiration (kg C m day
).
Heterotrophic respiration (kg C m day
).
Non-structural carbon (kg C m).
Biomass of seeds (kg C m).
Biomass of leaves (kg C m).
Biomass of fine roots (kg C m).
Biomass of sapwood (kg C m).
biomass of heartwood (kg C m).
Non-structural N pool (kg N m).
Nitrogen of seeds (kg N m).
Nitrogen of leaves (kg N m).
Nitrogen of roots (kg N m).
Nitrogen of sapwood (kg N m).
Nitrogen of heartwood (kg N m).
Microbial carbon (kg C m).
Fast soil carbon pool (kg C m).
Slow soil carbon pool (kg C m).
Microbial nitrogen (kg N m).
Fast soil nitrogen pool (kg N m).
Slow soil nitrogen pool (kg N m).
Mineral nitrogen pool (kg N m).
Nitrogen uptake (kg N m).
output_daily_cohorts
A data.frame with daily predictions for each canopy cohort.
Year of the simulation.
Day of the year.
Hour of the day.
An integer indicating the cohort identity.
An integer indicating the Plant Functional Type.
An integer indicating the crown layer, numbered from top to bottom.
Number of trees per area (trees ha).
Fraction of layer area occupied by this cohort.
Leaf area index (m/m
).
Gross primary productivity (kg C tree day
).
Plant autotrophic respiration (kg C tree day
).
Transpiration (mm tree day
).
Carbon allocated to leaves (kg C tree day
).
Carbon allocated to fine roots (kg C tree day
).
Carbon allocated to wood (kg C tree day
).
Nonstructural carbohydrates of a tree in this cohort (kg C
tree).
Seed biomass of a tree in this cohort (kg C tree).
Leaf biomass of a tree in this cohort (kg C tree).
Fine root biomass of a tree in this cohort (kg C tree).
Sapwood biomass of a tree in this cohort (kg C tree).
Heartwood biomass of a tree in this cohort (kg C tree).
Nonstructural nitrogen of a tree in this cohort (kg N tree).
Seed nitrogen of a tree in this cohort (kg N tree).
Leaf nitrogen of a tree in this cohort (kg N tree).
Fine root nitrogen of a tree in this cohort (kg N tree).
Sapwood nitrogen of a tree in this cohort (kg N tree).
Heartwood nitrogen of a tree in this cohort (kg N tree).
output_annual_tile
A data.frame with annual outputs at tile level.
Year of the simulation.
Crown area index (m/m
).
Leaf area index (m/m
).
Number of trees per area (trees ha).
Diameter at tile level (cm).
Tree density for trees with DBH > 12 cm (individuals
ha).
Diameter at tile level considering trees with DBH > 12 cm (cm).
Quadratic mean diameter at tile level considering trees with DBH > 12 cm (cm).
Net primary productivity (kg C m yr
).
Gross primary productivity (kg C m yr
).
Plant autotrophic respiration (kg C m yr
).
Heterotrophic respiration (kg C m yr
).
Annual precipitation (mm m yr
).
Soil water content in root zone (kg m).
Transpiration (mm m yr
).
Evaporation (mm m yr
).
Water runoff (mm m yr
).
Plant biomass (kg C m).
Soil carbon (kg C m).
Plant nitrogen (kg N m).
Soil nitrogen (kg N m).
Total nitrogen in plant and soil (kg N m).
Nonstructural carbohydrates (kg C m).
Seed biomass (kg C m).
Leaf biomass (kg C m).
Fine root biomass (kg C m).
Sapwood biomass (kg C m).
Heartwood biomass (kg C m).
Nonstructural nitrogen (kg N m).
Seed nitrogen (kg N m).
Leaf nitrogen (kg N m).
Fine root nitrogen (kg N m).
Sapwood nitrogen (kg N m).
Heartwood nitrogen (kg N m).
Microbial carbon (kg C m).
Fast soil carbon pool (kg C m).
Slow soil carbon pool (kg C m).
Microbial nitrogen (kg N m).
Fast soil nitrogen pool (kg N m).
Slow soil nitrogen pool (kg N m).
Mineral nitrogen pool (kg N m).
Nitrogen fixation (kg N m).
Nitrogen uptake (kg N m).
Annual available nitrogen (kg N m).
Annual nitrogen from plants to soil (kg N m).
Annual nitrogen loss (kg N m).
Total seed carbon (kg C m).
Total seed nitrogen (kg N m).
Total carbon from all compartments but seeds
(kg C m).
Total nitrogen from all compartments but seeds
(kg N m).
Age of the oldest tree in the tile (years).
Maximum volumne of a tree in the tile (m).
Maximum DBH of a tree in the tile (m).
Growth of a tree, including carbon allocated to leaves
(kg C m year
).
Growth of a tree, including carbon allocated to sapwood
(kg C m year
).
Number of trees that died (trees m year
).
Carbon biomass of trees that died (kg C
m year
).
Continuous biomass turnover (kg C m year
).
Carbon turnover rate, calculated as the ratio
between plant biomass and NPP (year).
output_annual_cohorts
A data.frame of annual outputs at the cohort level.
Year of the simulation.
An integer indicating the cohort identity.
An integer indicating the Plant Functional Type.
An integer indicating the crown layer, numbered from top to bottom.
Number of trees per area (trees ha).
Fraction of layer area occupied by this cohort.
Diameter growth of a tree in this cohort (cm year).
Tree diameter (cm).
Tree height (m).
Age of the cohort (years).
Crown area of a tree in this cohort (m).
Sum of sapwood and heartwood biomass of a tree in this cohort
(kg C tree).
Nonstructural carbohydrates in a tree (kg C tree).
Nonstructural nitrogen of a tree (kg N tree).
Total growth of a tree, including carbon allocated to seeds,
leaves, fine roots, and sapwood (kg C tree year
).
Fraction of carbon allocated to seeds to total growth.
Fraction of carbon allocated to leaves to total growth.
Fraction of carbon allocated to fine roots to total growth.
Fraction of carbon allocated to sapwood to total growth.
Gross primary productivity of a tree (kg C tree
year
).
Net primary productivity of a tree (kg C tree
year
).
Plant autotrophic respiration (kg C tree yr
).
Nitrogen uptake (kg N tree yr
).
Nitrogen fixation (kg N tree yr
).
Maximum leaf area index for a tree (m m
).
Tree volume (m).
Number of trees that died (trees yr).
Carbon biomass of trees that died (kg C yr).
Mortality rate of this cohort (yr).
# Example BiomeE model run # Use example drivers data drivers <- biomee_gs_leuning_drivers # Run BiomeE for the first site mod_output <- run_biomee_f_bysite( sitename = drivers$sitename[1], params_siml = drivers$params_siml[[1]], site_info = drivers$site_info[[1]], forcing = drivers$forcing[[1]], params_tile = drivers$params_tile[[1]], params_species = drivers$params_species[[1]], init_cohort = drivers$init_cohort[[1]], init_soil = drivers$init_soil[[1]] )
# Example BiomeE model run # Use example drivers data drivers <- biomee_gs_leuning_drivers # Run BiomeE for the first site mod_output <- run_biomee_f_bysite( sitename = drivers$sitename[1], params_siml = drivers$params_siml[[1]], site_info = drivers$site_info[[1]], forcing = drivers$forcing[[1]], params_tile = drivers$params_tile[[1]], params_species = drivers$params_species[[1]], init_cohort = drivers$init_cohort[[1]], init_soil = drivers$init_soil[[1]] )
Run P-model Fortran model on single site.
run_pmodel_f_bysite( sitename, params_siml, site_info, forcing, params_modl, makecheck = TRUE, verbose = TRUE )
run_pmodel_f_bysite( sitename, params_siml, site_info, forcing, params_modl, makecheck = TRUE, verbose = TRUE )
sitename |
Site name. |
params_siml |
Simulation parameters. |
site_info |
Site meta info in a data.frame. |
forcing |
A data frame of forcing climate data, used as input. |
params_modl |
A named list of free (calibratable) model parameters. See |
makecheck |
A logical specifying whether checks are performed
to verify forcings and model parameters. |
verbose |
A logical specifying whether to print warnings.
Defaults to For further specifications of above inputs and examples see |
Depending on the input model parameters, it's possible to run the different P-model setups presented in Stocker et al. 2020 GMD. The P-model version implemented in this package allows more flexibility than the one presented in the paper, with the following functions:
The temperature dependence of the quantum yield efficiency is given by: if
,
if
, and
if
.
The ORG setup can be reproduced by setting kphio_par_a = 0
and calibrating the kphio
parameter only.
The BRC setup (which calibrates in Eq. 18) is more difficult to reproduce,
since the temperature-dependency has been reformulated and a custom cost
function would be necessary for calibration. The new parameters
are related to
as follows:
The soil moisture stress is implemented as if
and
if
.
In Stocker et al. 2020 GMD, the threshold plant-available soil water is set as
= 0.6 * whc
where whc
is the site's water holding capacity. Also,
the reduction at low soil moisture (
) was parameterized
as a linear function of mean aridity (Eq. 20 in Stocker et al. 2020 GMD) but is
considered a constant model parameter in this package.
Hence, the FULL calibration setup cannot be
exactly replicated.
Model output is provided as a tidy dataframe, with columns:
date
Date of the observation in YYYY-MM-DD format.
year_dec
Decimal representation of year and day of the year (for example, 2007.000 corresponds to 2007-01-01 and 2007.003 to 2007-01-02.
fapar
Fraction of photosynthetic active radiation (fAPAR), taking values between 0 and 1.
gpp
Gross Primary Productivity (GPP) for each time stamp
(in gC m d
).
aet
Actual evapotranspiration (AET), calculated by SPLASH following Priestly-Taylor (in mm d).
le
Latent heat flux (in J m d
).
pet
Potential evapotranspiration (PET), calculated by SPLASH following Priestly-Taylor (in mm d).
vcmax
Maximum rate of RuBisCO carboxylation
(Vcmax) (in mol C m d
).
jmax
Maximum rate of electron transport for RuBP regeneration
(in mol CO m
s
).
vcmax25
Maximum rate of carboxylation (Vcmax),
normalised to 25C (in mol C m
d
).
jmax25
Maximum rate of electron transport, normalised to
25C (in mol C m
s
).
gs_accl
Acclimated stomatal conductance (in
mol C m d
Pa
).
wscal
Relative soil water content, between 0 (permanent wilting point, PWP) and 1 (field capacity, FC).
chi
Ratio of leaf-internal to ambient CO, ci:ca (unitless).
iwue
Intrinsic water use efficiency (iWUE) (in Pa).
rd
Dark respiration (Rd) in gC m d
.
tsoil
Soil temperature, in C.
netrad
Net radiation, in W m. WARNING: this is currently ignored as a model forcing. Instead, net radiation is internally calculated by SPLASH.
wcont
Soil water content, in mm.
snow
Snow water equivalents, in mm.
cond
Water input by condensation, in mm d
# Define model parameter values from previous work params_modl <- list( kphio = 0.04998, # setup ORG in Stocker et al. 2020 GMD kphio_par_a = 0.0, # disable temperature-dependence of kphio kphio_par_b = 1.0, soilm_thetastar = 0.6 * 240, # old setup with soil moisture stress soilm_betao = 0.0, beta_unitcostratio = 146.0, rd_to_vcmax = 0.014, # from Atkin et al. 2015 for C3 herbaceous tau_acclim = 30.0, kc_jmax = 0.41 ) # Run the Fortran P-model mod_output <- run_pmodel_f_bysite( # unnest drivers example data sitename = p_model_drivers$sitename[1], params_siml = p_model_drivers$params_siml[[1]], site_info = p_model_drivers$site_info[[1]], forcing = p_model_drivers$forcing[[1]], params_modl = params_modl )
# Define model parameter values from previous work params_modl <- list( kphio = 0.04998, # setup ORG in Stocker et al. 2020 GMD kphio_par_a = 0.0, # disable temperature-dependence of kphio kphio_par_b = 1.0, soilm_thetastar = 0.6 * 240, # old setup with soil moisture stress soilm_betao = 0.0, beta_unitcostratio = 146.0, rd_to_vcmax = 0.014, # from Atkin et al. 2015 for C3 herbaceous tau_acclim = 30.0, kc_jmax = 0.41 ) # Run the Fortran P-model mod_output <- run_pmodel_f_bysite( # unnest drivers example data sitename = p_model_drivers$sitename[1], params_siml = p_model_drivers$params_siml[[1]], site_info = p_model_drivers$site_info[[1]], forcing = p_model_drivers$forcing[[1]], params_modl = params_modl )
Runs BiomeE model for multiple sites.
runread_biomee_f(drivers, makecheck = TRUE, parallel = FALSE, ncores = 2)
runread_biomee_f(drivers, makecheck = TRUE, parallel = FALSE, ncores = 2)
drivers |
A nested data frame with one row for each site and columns
named according to the arguments of function |
makecheck |
A logical specifying whether checks are performed
to verify forcings and model parameters. |
parallel |
Flag specifying whether simulations are to be
parallelised (sending data from a certain number of sites to each core).
Defaults to |
ncores |
An integer specifying the number of cores used for parallel computing. Defaults to 2. |
A data frame (tibble) with one row for each site, site information
stored in the nested column site_info
and model outputs stored in the
nested column data
. See run_biomee_f_bysite
for a detailed
description of the outputs.
Example outputs are provided as p_model_output
and
p_model_output_vcmax25
.
# Example BiomeE model run runread_biomee_f( drivers = biomee_gs_leuning_drivers ) runread_biomee_f( drivers = biomee_p_model_drivers )
# Example BiomeE model run runread_biomee_f( drivers = biomee_gs_leuning_drivers ) runread_biomee_f( drivers = biomee_p_model_drivers )
Runs P-model for multiple sites.
runread_pmodel_f(drivers, par, makecheck = TRUE, parallel = FALSE, ncores = 1)
runread_pmodel_f(drivers, par, makecheck = TRUE, parallel = FALSE, ncores = 1)
drivers |
A nested data frame with one row for each site and columns
named according to the arguments of function |
par |
A named list of free (calibratable) model parameters.
|
makecheck |
A logical specifying whether checks are performed
to verify forcings and model parameters. |
parallel |
A logical specifying whether simulations are to be
parallelised (sending data from a certain number of sites to each core).
Defaults to |
ncores |
An integer specifying the number of cores used for parallel
computing (by default |
Depending on the input model parameters, it's possible to run the different P-model setups presented in Stocker et al. 2020 GMD. The P-model version implemented in this package allows more flexibility than the one presented in the paper, with the following functions:
The temperature dependence of the quantum yield efficiency is given by: if
,
if
, and
if
.
The ORG setup can be reproduced by setting kphio_par_a = 0
and calibrating the kphio
parameter only.
The BRC setup (which calibrates in Eq. 18) is more difficult to reproduce,
since the temperature-dependency has been reformulated and a custom cost
function would be necessary for calibration. The new parameters
are related to
as follows:
The soil moisture stress is implemented as if
and
if
.
In Stocker et al. 2020 GMD, the threshold plant-available soil water is set as
= 0.6 * whc
where whc
is the site's water holding capacity. Also,
the reduction at low soil moisture (
) was parameterized
as a linear function of mean aridity (Eq. 20 in Stocker et al. 2020 GMD) but is
considered a constant model parameter in this package.
Hence, the FULL calibration setup cannot be
exactly replicated.
A data frame (tibble) with one row for each site, site information
stored in the nested column site_info
and outputs stored in the nested
column data
. See run_pmodel_f_bysite
for a detailed
description of the outputs.
Example outputs are provided as biomee_p_model_output
and
biomee_gs_leuning_output
.
# Define model parameter values from previous work params_modl <- list( kphio = 0.04998, # setup ORG in Stocker et al. 2020 GMD kphio_par_a = 0.0, # disable temperature-dependence of kphio kphio_par_b = 1.0, soilm_thetastar = 0.6 * 240, # old setup with soil moisture stress soilm_betao = 0.0, beta_unitcostratio = 146.0, rd_to_vcmax = 0.014, # from Atkin et al. 2015 for C3 herbaceous tau_acclim = 30.0, kc_jmax = 0.41 ) # Run the model for these parameters and the example drivers output <- rsofun::runread_pmodel_f( drivers = rsofun::p_model_drivers, par = params_modl) output_vcmax25 <- rsofun::runread_pmodel_f( drivers = rsofun::p_model_drivers_vcmax25, par = params_modl)
# Define model parameter values from previous work params_modl <- list( kphio = 0.04998, # setup ORG in Stocker et al. 2020 GMD kphio_par_a = 0.0, # disable temperature-dependence of kphio kphio_par_b = 1.0, soilm_thetastar = 0.6 * 240, # old setup with soil moisture stress soilm_betao = 0.0, beta_unitcostratio = 146.0, rd_to_vcmax = 0.014, # from Atkin et al. 2015 for C3 herbaceous tau_acclim = 30.0, kc_jmax = 0.41 ) # Run the model for these parameters and the example drivers output <- rsofun::runread_pmodel_f( drivers = rsofun::p_model_drivers, par = params_modl) output_vcmax25 <- rsofun::runread_pmodel_f( drivers = rsofun::p_model_drivers_vcmax25, par = params_modl)