In miquelcaceres/IFNdyn: IFN forest dynamics and scenario analysis

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

Introduction

Package IFNdyn was designed as a tool to model forest dynamics and management effects on permanent forest plots of the Spanish National Forest Inventory (Inventario Forestal Nacional). The package implements tree-level distance-independent growth and yield models, including growth, survival and recruitment processes. The package includes functions to simulate forest dynamics and silvicultural treatments, and ancillary functions to estimate wood volume and biomass for different plant components. Importantly, through its high-level function called IFNscenario(), the package allows conducting large-scale simulations of forest dynamics under user-defined management and climate scenarios. This document introduces the design of the different models of forest dynamics and provides a user guide for the functions included in the package.

library(IFNdyn)
library(knitr)
library(ggplot2)
data("species_selection")

Input data

Tree records

The basic data structure of the IFNdyn package is a data frame with as many rows as tree records and columns:

• ID: Plot code (unique within IFN), where the first one/two digit(s) indicate the province and the remaining indicate the plot number within the province.
• Species: Tree species code.
• N: Tree density (in number of individuals per hectare) that the record represents.
• DBH: Tree diameter at breast height (cm).
• H: Tree height (m).

The data set exampleTreeData is provided as an example of such structure (it includes additional columns not relevant for users of the package):

data("exampleTreeData")
head(exampleTreeData)

Plot data

Most package functions require input data of variables measured at forest plot level. Specifically, a data frame is often required with forest plot in rows (row names are unique forest plot codes) and columns:

• Province: Spanish province code (1 to 50).
• Rad: Annual mean of daily solar radiation (MJ/day).
• Temp: Annual mean of daily mean temperature (degrees C).
• Prec: Mean annual precipitation (mm).
• PET: Mean annual potential evapo-transpiration (mm).
• elevation: Elevation above sea level (m).
• slope: Slope in degrees.
• SWHC: Soil water holding capacity (mm).
• Managed: Boolean to indicate whether the plot can be managed or not. These can be used, for example, to avoid simulating silviculture in forest plots that are within the bounds of protected areas (only needed for IFNscenario()).
• Estrato: Stratum code. IFN strata identify areas covered by more or less homogeneous forests. They are like forest classes used extrapolate the results of calculations done on forest plots to larger areas (only needed for IFNscenario()).

As many users may find difficult to assemble environmental data, the package includes a data set called plotDataIFN with precalculated values for IFN3 plots:

data("plotDataIFN")
head(plotDataIFN)

Spatial coordinates are needed to model ingrowth including species colonization. The package includes the coordinates for three example forest plots only:

data(examplePlotCoords)
examplePlotCoords

Models of forest dynamics

The IFNdyn package a number of functions to simulate population processes (growth, survival and ingrowth) individually as well as a function to simulate all processes sequentially (i.e. forest dynamics). The functions implement empirical models that have been calibrated using data from the second, third and fourth Spanish Forest Inventories (IFN2, IFN3, IFN4). Since Spanish forest inventory surveys are typically repeated at 10-yr intervals, the models are best suited for simulations using this time step, but they can be used for smaller (e.g. annual) time steps.

When developing IFN models, tree species occurring in Peninsular Spain, Balearic Islands and Canary Islands were grouped to decrease the number of models to be calibrated (otherwise, there would not be enough records for all species): r nrow(species_selection) species and species groups were defined for model development. The following table relates the species codes with the species/group name:

df = species_selection[,c(2,1)]
df$names = as.character(df$names)
df$species = as.character(df$species)
names(df) <-c("Species/group name", "Species codes")
df[nrow(df),1] = paste0(substr(df[nrow(df),1],1,60),"...")
kable(df, format = "html", digits=2)

Note that trying to run forest dynamics on plot records with species codes not included in this list will result in a execution error. Six kinds of models are defined for each species or species group (available in functions given in brackets):

1. Individual tree section increment (cm$^2$) or diameter increment (cm) [IFNgrowth()].
2. Individual tree heigh increment (m) [IFNheightgrowth()].
3. Survival probability (%) of an individual tree [IFNsurvival() and IFNsurvivalPG()].
4. Ingrowth density corresponding to the diameter class between 7.5 and 12.4 cm [IFNingrowth() and IFNingrowthdisp()].
5. Average diameter (cm) of trees entering the diameter class between 7.5 and 12.4 cm [IFNingrowth() and IFNingrowthdisp()].
6. Total height of an individual tree (m) [IFNheight()].

Explanatory variables

Models of forest dynamics use a set of predictors defined at the tree or stand levels [R representation in brackets]:

a. Tree diameter at breast height (cm) [d] b. Inverse of tree diameter at breast height (cm$^{-1}$) [inv.d] c. Logarithm of tree diameter at breast height [ln.d] d. Square of tree diameter at breast height (cm$^{2}$) [sq.d] e. Square-root of tree diameter at breast height (cm$^{1/2}$) [sqrt.d] f. Tree height (m) [d] g. Inverse of tree height (m$^{-1}$) [inv.h] h. Logarithm of tree height [ln.h] i. Basal area of trees larger than the target tree (m$^2$·ha$^{-1}$) [BAL] j. Basal area of larger trees divided by the logarithm of tree diameter [BAL.ln.d] k. Basal area of trees extracted in a previous step (not normally used; m$^2$·ha$^{-1}$) [BAL.ext] l. Stand total basal area (m$^2$·ha$^{-1}$) [G] m. Logarithm of total basal area [ln.G] n. Logarithm of stand tree density [ln.N] o. Proportion of trees of the target species [Nsp.N] p. Standard deviation of diameter at breast height [sd.DBH]

The predictors listed above are calculated from a data frame of tree records using functions prepareTreeCompetitionTable(), prepareSpeciesCompetitionTable() and prepareSpatialSpeciesCompetitionTable(). Users of the package will not normally use these three functions directly, but we illustrate them here because it is important to understand that they are used internally as intermediate steps in forest dynamics simulations.

Diameter/height increment and survival models require the following tree and stand explanatory variables:

head(prepareTreeCompetitionTable(exampleTreeData, provinceFromID = TRUE))

Ingrowth models (local recruitment) require somewhat different variables, particularly Nsp.N:

prepareSpeciesCompetitionTable(exampleTreeData, provinceFromID = TRUE)

Additional information is required as input for new species ingrowth (i.e. species colonization). First, one needs an object of neighbours and their distances can be build from a matrix of spatial coordinates using function IFNknn() (here restricted to one neighbour):

nb = IFNknn(examplePlotCoords, k = 1, verbose = FALSE)
nb

To prepare input for colonization we also need to know the relative dominance of the target species [Nsp.N] in all possible neighbour plots (possibly including plots different from the target plots). Here we build such a matrix corresponding to the three example plots:

Nmatrix = tapply(exampleTreeData$N, list(exampleTreeData$ID, exampleTreeData$Species),
                 sum, na.rm=T)
NspNmatrix = sweep(Nmatrix, 1, rowSums(Nmatrix, na.rm=T), "/")
NspNmatrix = rbind(NspNmatrix,c(0.5, NA, NA, NA))
rownames(NspNmatrix)[4] = "80003"
NspNmatrix = cbind(NspNmatrix, c(NA,NA,NA,0.5))
colnames(NspNmatrix)[5] = "25"
NspNmatrix

The function to prepare data input for ingrowth models including colonization effects is prepareSpatialSpeciesCompetitionTable():

prepareSpatialSpeciesCompetitionTable(exampleTreeData, 
                                      nb = nb, NspNmatrix = NspNmatrix,
                                      provinceFromID = TRUE)

Note that the first six rows [where Nsp.N > 0] are the same as those obtained with prepareSpeciesCompetitionTable(), referring to local recruitment. The last three rows [where Nsp.N == 0] are used to determine colonizations, which require additional explanatory variables [Nsp.N_s0.5 to Nsp.N_s4].

Function drawSubmodelEffects() allows examining the effect of different predictors on the response variable of the different submodels. This function will be illustrated below, when describing each of the submodels.

Diameter increment

Tree diameter increment is estimated using a logarithmic function with multiplicative effects:

$$\ln(DI) = a_0 + b_1 \cdot d+b_2/d + b_3 \cdot \ln(d+1) + b_4 \cdot d^2 + b_5 \cdot \sqrt{d}\ +b_6 \cdot \ln(G) + b_7 \cdot BAL+b_8 \cdot (BAL/\ln(d+1))+b_{9} \cdot BAL_{ext}\ +b_{10} \cdot (Temp-15)^2 + b_{11} \cdot (P/PET) + b_{12} \cdot SWHC + b_{13} \cdot Rad + b_{14} \cdot slope$$

Parameters for each species were estimated by fitting Generalized Linear Models with Gamma response and a logarithmic function as link. They are stored in a data frame called defaultGrowthParams:

data("defaultGrowthParams")
head(defaultGrowthParams[,1:17])

The following code allows showing the effect of initial DBH on tree diameter increment for six different pine species:

pines = c("21", "22", "23", "24", "25", "26")
drawSubmodelEffects(species = pines, predictor = "DBH", submodel = "growth")

Analogous plots can be drawn to inspect the effect of stand basal area ($G$), or the basal area of larger trees ($BAL$):

drawSubmodelEffects(species = pines, predictor = "G", submodel = "growth")

drawSubmodelEffects(species = pines, predictor = "BAL", submodel = "growth")

Function IFNgrowth() allows calculating diameter increments for an input data frame with tree records, as long as plot environmental conditions are also provided:

IFNgrowth(exampleTreeData, plotDataIFN)

By default, the function returns increments expected for a 10-yr period, but other growth intervals (in years) can be specified.

Height increment

Height diameter increment is also estimated using a logarithmic function with multiplicative effects:

$$\ln(HI) = a_0 + b_2/h + b_2 \cdot \ln(h+1) + b_3 \cdot (h/d)\ + b_{4} \cdot G +b_5 \cdot (BAL/\ln(d+1))+b_{6} \cdot BAL_{ext} \ +b_{7} \cdot (Temp-15)^2 + b_{8} \cdot (P/PET) + b_{9} \cdot SWHC + b_{10} \cdot Rad + b_{11} \cdot slope$$

Parameters for each species were estimated by fitting Generalized Linear Models with Gamma response and a logarithmic function as link. They are stored in a data frame called defaultHeightGrowthParams:

data("defaultHeightGrowthParams")
head(defaultHeightGrowthParams[,1:15])

The following code allows showing the effect of initial DBH on tree diameter increment for six different pine species:

pines = c("21", "22", "23", "24", "25", "26")
drawSubmodelEffects(species = pines, predictor = "H", submodel = "heightgrowth")

Analogous plots can be drawn to inspect the effect of stand basal area ($G$):

drawSubmodelEffects(species = pines, predictor = "G", submodel = "heightgrowth")

Analogously to diameter increment, function IFNheightgrowth() allows calculating diameter increments for an input data frame with tree records:

IFNheightgrowth(exampleTreeData, plotDataIFN)

By default, the function returns increments expected for a 10-yr period, but other growth intervals (in years) can be specified.

Static height

The height increment model cannot always be applied (it is not available for some species and cannot be used to estimate the height of ingrowth). In this case, heights are updated according to a static height submodel. The following non-linear model form is used:

$$H = \frac{a_0+b_1 \cdot (P/PET) +b_2 \cdot SWHC +b_3 \cdot Rad + b_5 \cdot elevation}{1+(c_1/d)+(c_2/d^2)}$$

Parameters for each species were estimated by fitting non-linear models. They are stored in a data frame called defaultHeightParams:

data("defaultHeightParams")
head(defaultHeightParams[,1:9])

The following figures show the effect of DBH and (P/PET) on tree height:

drawSubmodelEffects(species = pines, predictor = "DBH", submodel = "height")

drawSubmodelEffects(species = pines, predictor = "P/PET", submodel = "height")

Function IFNheight() returns the height (in m) for an input data frame with tree records:

IFNheight(exampleTreeData, plotDataIFN)

Survival

Probability of tree survival is modelled using the following logistic model:

$$logit(P_{surv}) = a_0 + b_1 \cdot d+b_2 \cdot d^2+b_3 \cdot \sqrt{d}+b_4 \cdot h \ +b_5 \cdot BAL+b_6 \cdot BAL_{ext} +b_7 \cdot \ln(G)\ +b_{8} \cdot (Temp-15)^2 + b_{9} \cdot (P/PET) + b_{10} \cdot SWHC + b_{11} \cdot Rad + b_{12} \cdot slope$$

Parameters for each species were estimated by fitting Generalized Linear Models with binomial response and a logit link function. They are stored in a data frame called defaultSurvivalParams:

data("defaultSurvivalParams")
head(defaultSurvivalParams)

The following figures show the effect of tree DBH and the basal area of larger trees on survival probability for pine species:

drawSubmodelEffects(species = pines, predictor = "DBH", submodel = "survival")

drawSubmodelEffects(species = pines, predictor = "BAL", submodel = "survival")

Function IFNsurvival() allows calculating survival probability for an input data frame with tree records:

IFNsurvival(exampleTreeData, plotDataIFN)

By default, the function returns probabilities expected for a 10-yr period, but other time steps (in years) can be specified. An alternative function is available, IFNsurvivalPG() that allows including the diameter increment in a previous time step as predictor of survival.

Ingrowth (local recruitment)

Ingrowth is modelled as the incorporation of new trees to the diametric class [7.5-12.5] cm. The following equation is used to estimate the presence of ingrowth:

$$logit(B_{new}) = a_0 + b_1 \cdot \ln(N)+b_2 \cdot (N_{sp}/N)+b_3 \cdot \ln{G}+b_4 \cdot SD(d)\ +b_5 \cdot (Temp-15)^2 + b_6 \cdot (P/PET) + b_7 \cdot SWHC + b_8 \cdot Rad + b_9 \cdot slope$$

The following equation is used to estimate the number (density) of new trees:

$$\ln(N_{new}) = a_0 + b_1 \cdot \ln(N)+b_2 \cdot (N_{sp}/N)+b_3 \cdot \ln{G}+b_4 \cdot SD(d)\ +b_5 \cdot (Temp-15)^2 + b_6 \cdot (P/PET) + b_7 \cdot SWHC + b_8 \cdot Rad + b_9 \cdot slope$$ Parameters for each species were estimated by fitting Generalized Linear Models with Binomial and Gamma response and a logarithmic function as link, respectively. They are stored in a data frames called defaultIncorporationBParams and defaultIncorporationNParams:

data("defaultIngrowthBParams")
head(defaultIngrowthBParams)
data("defaultIngrowthNParams")
head(defaultIngrowthNParams)

Presence and density of ingrowth are calculated simultaneously. The following figures illustrate the effect of stand basal area and density on the final density of ingrowth:

drawSubmodelEffects(species = pines, predictor = "G", submodel = "ingrowth")

drawSubmodelEffects(species = pines, predictor = "N", submodel = "ingrowth")

Tree diameter of new trees within the diametric class [7.5-12.5] is modelled using:

$$\ln(DBH_{new}) = a_0 + b_1 \cdot G + b_2 \cdot CSWD + b_3 \cdot Rad + b_4 \cdot Temp + b_5 \cdot slope$$

Parameters for each species were estimated by fitting Generalized Linear Models with Gamma response and a logarithmic function as link. They are stored in a data frame called defaultIncorporationDBHParams:

data("defaultIngrowthDBHParams")
head(defaultIngrowthDBHParams)

Function IFNingrowth() returns a data frame with ingrowth records (including density, DBH and height) for an input data frame with tree records:

IFNingrowth(exampleTreeData, plotDataIFN)

Like growth and survival functions, by default IFNingrowth() determines the ingrowth expected for a 10-yr period, but other time steps (in years) can be specified.

Ingrowth (colonization)

Ingrowth presence, density and diameter corresponding to species colonization are modelled in the same way as for local recruitment, but replacing the role of [Nsp.N] by a set of variables with different weight of neighbours [Nsp.N_s0.5 to Nsp.N_s4]. Parameters are stored in a data frames called defaultIncorporationBParams_disp, defaultIncorporationNParams_disp and defaultIncorporationDBHParams_disp. As explained above, determining colonization requires identifying the neighbours of each plot and the relative abundance of the target species in those neighbours (i.e. objects nb and Nsp.Nmatrix created above). Once all the necessary data is available, function IFNingrowthdisp() returns a data frame of new species ingrowth (with density, DBH and height) for an input set of forest stands:

IFNingrowthdisp(exampleTreeData, plotDataIFN, 
                nb = nb, NspNmatrix)

As before, by default IFNingrowthdisp() determines the new species colonization expected for a 10-yr period.

Running all submodels together

Function IFNdynamics() allows conducting a simulation of forest dynamics for a number of time steps (by default, each time step is 10 years and one time step is simulated).

treeDataNew = IFNdynamics(exampleTreeData, plotDataIFN)

The following code shows the tree data frame of a single plot before and after simulation:

exampleTreeData[exampleTreeData$ID == "80001",]
treeDataNew[treeDataNew$ID == "80001",]

Note that the result of calling IFNdynamics() has different columns than the original data frame, which can contain extra columns not used in simulations. By default, species names are added in the output, for convenience (see function speciesNames()).

Function IFNdynamics() has several options. In particular, we can simulate dynamics for 20 steps (i.e. 200 years) and keep the transient states, the latter by specifying sequence = TRUE:

res = IFNdynamics(exampleTreeData, plotDataIFN, numSteps = 20, sequence = TRUE)

Function standDynamicsSummary() allows drawing the result of long-term simulations. For example, below we draw the variation of stand basal area during the previous 200-yr simulation:

standDynamicsSummary(res, standVariable = "BA")

Another option is to run the models in stochastic mode (instead of deterministic mode)

res = IFNdynamics(exampleTreeData, plotDataIFN, numSteps = 20, sequence = TRUE, stochastic = TRUE)

which in our example produces the following output:

standDynamicsSummary(res, standVariable = "BA")

By default, IFNdynamics() uses the static height model, but we can force using height increment models setting the logical parameter heightgrowth to TRUE. Analogously, default runs of IFNdynamics() include local recruitment only, but we can use parameters nb and NspNmatrix to force simulation of species colonization. Other options can be inspected in the help page of the function.

Forest management

Package IFNdyn allows simulating two main kinds of silvicultural models:

1. Irregular: Include stand thinnings every time basal area exceeds a given threshold or the Hart-Becking index goes below a given threshold.
2. Regular: Includes both stand thinnings as well as final cuts and planting of new trees (saplings).

Silvicultural prescriptions

Silvicultural prescriptions need to be specified for each forest stand. Thus, the data structure for prescription is a data frame with forest plots in rows (row names are plot IDs) and a set of columns with prescription specifications. Prescription for thinnings require the basal area or Hart-Becking index thresholds (columns thinningBA and thinningHB), the thinning intensity (i.e. the percentage of basal area extracted, column thinningPerc) and the thinning type (column thinning), which can be of the following kinds:

• above: The largest trees are removed until the target basal area to be extracted is reached.
• below: The smallest trees are removed until the target basal area to be extracted is reached.
• intermediate: Trees from all diameter classes are extracted with equal probability until the target basal area to be extracted is reached.
• below-intermediate: Half of the basal area to be extracted is cut as intermediate and the other half as below.
• above-intermediate: Half of the basal area to be extracted is cut as intermediate and the other half as above.
• user-defined: The user provides diameter thresholds and percentage of thinnings applied to each diameter class (see example below).

Prescription for final cuts require specifying the threshold for plot mean DBH that leads to final cuts when exceeded (column finalMeanDBH). Final cuts can be scheduled in one or more time steps. The user needs to specify the percentage of basal area to be removed each time step (column finalPerc). Finally, planting after final cuts is also specified, including the species to be planted, the initial density, initial DBH and initial height (columns plantingSpecies, plantingDBH, plantingHeight and plantingDensity).

An example of prescription data frame, corresponding to the forest plots in exampleTreeData, is given in object examplePrescriptionData:

data(examplePrescriptionData)
examplePrescriptionData

Column finalPreviousStage is used to store the action done in the previous time step, among those included in the final cuts (0 - before final cut; 1 - first final cut; 2 - second final cut; etc.). This is necessary for final cuts in regular models, which follow an order of treatments that cannot be inferred from the stand structure itself.

Management function

Given an input tree data frame and prescription data frame, function IFNmanagement() determines for each plot whether management has to be applied and, if so, it performs the silvicultural treatment:

res = IFNmanagement(exampleTreeData, examplePrescriptionData)

The function returns a list that includes three data frames. The first two are data frames treeDataOut and treeDataCut, which contain records for non-cut trees and cut trees, respectively:

head(res$treeDataOut)
head(res$treeDataCut)

Additionally, the result includes the silvicultural prescription for the next time step (with an update of finalPreviousStage, if needed):

res$nextPrescription

Simulation of demand-based scenarios

The main purposes underlying the development of IFNdyn was the need to simulate climate and/or management scenarios of forest dynamics over large areas. More specifically, this is achieved by stating a wood volume demand for the target area (typically a Spanish province) and performing simulations of forest dynamics on the set of forest plots within the area, while performing silvicultural treatments to extract, as far as possible, the amount of wood necessary to fill the input demand. At the same time, climate scenarios can be specified by providing values of climatic predictor for each time step. All these calculations are done using function IFNscenario(), which returns for each time step simulated the biomass, volume and density of trees that are standing, dead or extracted. The following subsections describe the inputs and algorithmic procedure underlying function IFNscenario().

Simulation inputs

Main inputs of IFNscenario() are:

a. Data frame of tree records (for the whole province). b. Data frame of plot data (for the plots in the whole province). c. Prescription data by species. d. Forest strata data frame. e. Demand data frame. f. Data frame of dynamic (i.e. climate predictors) data (optional). g. An object specifying plot neighbourhoods (optional).

Tree records can be obtained by reading tree data from a data base (see functions readPiesMayoresIFNX() in package IFNread). As explained above, the package includes a data set called plotDataIFN with precalculated values for IFN3 plots. Here, we describe the remaining three input items:

Wood-extraction demand

Wood demand is specified using the m3 of wood to be extracted for each species. The package provides default tables (one per Spanish province) with the 2005-2015 average cubic meters extracted by species, obtained from the 'Anuario de Estadística Forestal':

data("defaultDemand")
names(defaultDemand)

The following shows the demand table for the Lleida province:

head(defaultDemand$Lleida)

Column Species includes the species codes that are included and AnnualDemand contains the average annual wood volume extracted in cubic meters for the province.

Area covered by strata

IFNscenario() needs to know how large is the area represented by each forest plot included in the simulation. This is obtained from the classification of plots into 'Estratos'. The package includes information of the classification of forests of each Spanish province into strata, particularly the area covered by each one. For example, in Lleida province (25):

data(estratosIFN3)
head(estratosIFN3$Lleida)

Importantly, the row names of strata data frame need to match the column Estrato in the plot data frame (i.e. plotDataIFN):

head(plotDataIFN[plotDataIFN$Province=="25",])

Silvicultural prescription by species

In function IFNscenario() silvicultural prescription is the same as in function IFNmanagement(), except for the fact that it is not practical to specify the silvicultural model for each forest plot. This is solved by specifying the most usual model for each species or species group:

data("defaultPrescription")
head(defaultPrescription)

IFNscenario() determines the dominant species (in basal area) of each forest plot, and will apply the corresponding silvicultural model. defaultPrescription is a data frame with default silvicultural models for all tree species, but the user may modify it to better represent silvicultural practices in the target area.

Simulation algorithm

Function IFNscenario() performs the following operations:

1. Data preparation:
   • Removal of tree records with missing data.
   • Removal of stands with missing data.
   • Determination of dominant species and area represented for each forest plot.
   • Determination of non-managed plots according the species included in the demand data frame.
2. For each time step: a. Simulation of forest dynamics for non-managed plots. For each time step to be simulated:
   • Perform forest dynamics (growth, survival and ingrowth).
   • Store the volume and biomass of dead trees and standing trees. b. Simulation of forest dynamics and silvicultural treatments in managed plots. For each species/group in the demand:
   • Select the plots dominated by the target species.
   • Set the demand (in m3 per 10 years) to be fulfilled.
   • Apply silvicultural treatment to plots of the regular model that are at intermediate stages of final cut. Store the volume and biomass corresponding to cut trees and reduce the demand accordingly.
   • While demand is not fulfilled keep on treating the remaining plots. Plots are processed in order, according to the difference between the Hart-Becking index of the stand and the threshold indicated in the prescription. Each time a plot is treated, store the volume and biomass corresponding to the cut trees.
   • When demand is exceeded (or if no more plots are available), store the difference between the demand and the volume actually cut to modulate the demand of the next step.
   • Simulate forest dynamics (growth, survival and ingrowth) on the forest plots. Ingrowth is not simulated in those plots whose management model is regular, because planting is expected to occur at the end of the cycle.
   • Store the volume and biomass of dead trees and standing trees.
3. Assemble simulation results

knitr::include_graphics('fig_flowchart.png')

Post processing

The package allows merging the results of simulations on different target areas by calling function IFNscenario() on each of them and then calling mergeScenarios() with the results. The merging function will create an output with the same structure as the output of IFNscenario() but containing the results of all simulations together.

The package allows performing calculations of carbon sequestration/emissions with function scenarioEmissions(), which internally calls functions scenarioProducts(), scenarioDeadWood() and scenarioLiveCarbonStock().

Ancillary functions

Package IFNdyn provides several ancillary functions that underlie or complement the simulation of forest dynamics. In what follows, we describe them grouped by category.

Forest structure

A set of functions are included to perform different stand structure calculations from tree records:

• Function basalArea() returns basal area of stands, either the total or disaggregated by species.
• Function dominantHeight() returns the dominant height of stands.
• Function dominantDiameter() returns the dominant diameter of stands.
• Function hartBeckingIndex() returns the Hart-Becking index of stands.

As an example, the following code returns the basal area of each plot in exampleTreeData:

basalArea(exampleTreeData)

Tree volume and biomass

Calculations of biomass and volumes for trees and their components are obtained using the following functions:

• Function IFNbiomass() uses in-built empirical models to estimate the dry weight biomass (kg·ha$^{-1}$) of different plant fractions for each tree record.
• Function IFNvolume() uses in-built empirical models to estimate the volume of wood (in m$^3$·ha$^{-1}$) corresponding to each tree record.
• Function IFNproducts() takes a data frame of tree cuts and a table specifying product destination by species as input and calculates biomass or volume estimates corresponding to products and slash generated for each forest plot.

head(IFNbiomass(exampleTreeData))

head(IFNvolume(exampleTreeData, provinceFromID = TRUE))

head(IFNproducts(exampleTreeData, defaultProductsCAT))

Species names

The following two functions are used to obtain species names from IFN species codes:

• Function speciesNames() returns species names (i.e. binomial)
• Function speciesNamesModels() returns the name of the species or species group used in model development (see table above).
• Function functionalGroups() returns the functional group (conifer or hardwood) corresponding to the species names.

For example, the following code returns the species names corresponding to rows of the example tree data frame:

speciesNames(exampleTreeData$Species)

Additional processes

Perturbations are important drivers of forest dynamics. While different perturbations exist (landslides, windthrows, drought, fire, insect attacks...), IFNdyn includes only the possibility to simulate fire impacts in Catalonia, using function IFNfire(). The function returns two data frames, one containing the probability of fire in a 10-yr interval and the other with the probability of survival for each tree, if fire occurs.

fi = IFNfire(exampleTreeData, plotDataIFN)
fi$stand
head(fi$tree)

Function IFNfire() is called from IFNdynamics() or IFNscenario(), if the user wants to include fire impacts in simulations of forest dynamics.

miquelcaceres/IFNdyn documentation built on Feb. 1, 2021, 10:55 a.m.