Forest dynamics

Miquel De Caceres

2021-06-17

About this vignette

This document describes how to run the forest dynamics model of medfate. This document is meant to teach users to run the simulation model within R. Details of the model design and formulation can be found at https://emf-creaf.github.io/medfatebook/index.html.

Preparing model inputs

Any forest dynamics model needs information on climate, vegetation and soils of the forest stand to be simulated. Moreover, since models in medfate differentiate between species, information on species-specific model parameters is also needed. In this subsection we explain the different steps to prepare the data needed to run function fordyn().

Preparing model inputs

Model inputs are explained in greater detail in vignette ‘Simulation inputs’. Here we only review the different steps required to run function fordyn().

Soil, vegetation, meteorology and species data

Soil information needs to be entered as a data frame with soil layers in rows and physical attributes in columns. Soil physical attributes can be initialized to default values, for a given number of layers, using function defaultSoilParams():

spar = defaultSoilParams(2)

The soil input for water balance simulation is actually a list of class soil that is created using a function with the same name:

examplesoil = soil(spar)

As explained in the package overview, models included in medfate were primarily designed to be ran on forest inventory plots. Here we use the example object provided with the package:

data(exampleforestMED)
exampleforestMED
## $ID
## [1] "1"
## 
## $patchsize
## [1] 10000
## 
## $treeData
##   Species   N   DBH Height Z50  Z95
## 1      54 168 37.55    800 750 3000
## 2      68 384 14.60    660 750 3000
## 
## $shrubData
##   Species Cover Height Z50  Z95
## 1      65  3.75     80 300 1500
## 
## $herbCover
## [1] 10
## 
## $herbHeight
## [1] 20
## 
## attr(,"class")
## [1] "forest" "list"

In the basic water balance, only mean temperature, precipitation and potential evapotranspiration is required, but radiation may also be necessary to simulate snow melt.

data(examplemeteo)
head(examplemeteo)
##            MeanTemperature MinTemperature MaxTemperature Precipitation
## 2001-01-01      3.57668969     -0.5934215       6.287950      4.869109
## 2001-01-02      1.83695972     -2.3662458       4.569737      2.498292
## 2001-01-03      0.09462563     -3.8541036       2.661951      0.000000
## 2001-01-04      1.13866156     -1.8744860       3.097705      5.796973
## 2001-01-05      4.70578690      0.3288287       7.551532      1.884401
## 2001-01-06      4.57036721      0.5461322       7.186784     13.359801
##            MeanRelativeHumidity MinRelativeHumidity MaxRelativeHumidity
## 2001-01-01             78.73709            65.15411           100.00000
## 2001-01-02             69.70800            57.43761            94.71780
## 2001-01-03             70.69610            58.77432            94.66823
## 2001-01-04             76.89156            66.84256            95.80950
## 2001-01-05             76.67424            62.97656           100.00000
## 2001-01-06             89.01940            74.25754           100.00000
##            Radiation WindSpeed WindDirection       PET
## 2001-01-01  12.89251  2.000000           172 1.3212770
## 2001-01-02  13.03079  7.662544           278 2.2185985
## 2001-01-03  16.90722  2.000000           141 1.8045176
## 2001-01-04  11.07275  2.000000           172 0.9200627
## 2001-01-05  13.45205  7.581347           321 2.2914449
## 2001-01-06  12.84841  6.570501           141 1.7255058

Finally, simulations in medfate require a data frame with species parameter values, which we load using defaults for Catalonia (NE Spain):

data("SpParamsMED")

Simulation control

Apart from data inputs, the behaviour of simulation models can be controlled using a set of global parameters. The default parameterization is obtained using function defaultControl():

control = defaultControl("Granier")

Here we will run simulations of forest dynamics using the basic water balance model (i.e. transpirationMode = "Granier"). The complexity of the soil water balance calculations can be changed by using "Sperry" as input to defaultControl(). However, when running fordyn() sub-daily output will never be stored (i.e. setting subdailyResults = TRUE is useless).

Executing the forest dynamics model

In this vignette we will fake a three-year weather input by repeating the example weather data frame three times, while reducing precipitation of the second year to the half.

meteo2001 = examplemeteo
meteo2002 = examplemeteo
meteo2002$Precipitation = meteo2002$Precipitation/2
meteo2003 = examplemeteo
row.names(meteo2002) = seq(as.Date("2002-01-01"), 
                           as.Date("2002-12-31"), by="day")
row.names(meteo2003) = seq(as.Date("2003-01-01"), 
                           as.Date("2003-12-31"), by="day")
meteo_01_03 = rbind(meteo2001, meteo2002, meteo2003)

Now we run the forest dynamics model using all inputs (note that no intermediate input object is needed, as in spwb() or growth()):

fd<-fordyn(exampleforestMED, examplesoil, SpParamsMED, meteo_01_03, control, 
           latitude = 41.82592, elevation = 100)
## Simulating forest dynamics for year 2001 (1/3)
##    (a) Growth/mortality
##    (b) Recruitment
##        Coldest month mean temp. (Celsius): 0.05   Moisture index: 0.36   FPAR (%): 39.8
##        Tree species with seed rain: 54,68 recruited: <none>
##    (c) Summaries
## Simulating forest dynamics for year 2002 (2/3)
##    (a) Growth/mortality
##    (b) Recruitment
##        Coldest month mean temp. (Celsius): 0.05   Moisture index: 0.18   FPAR (%): 39.8
##        Tree species with seed rain: 54,68 recruited: <none>
##    (c) Summaries
## Simulating forest dynamics for year 2003 (3/3)
##    (a) Growth/mortality
##    (b) Recruitment
##        Coldest month mean temp. (Celsius): 0.05   Moisture index: 0.36   FPAR (%): 37.4
##        Tree species with seed rain: 54,68 recruited: <none>
##    (c) Summaries

It is worth noting that, while fordyn() calls function growth() internally for each simulated year, the verbose option of the control parameters only affects function fordyn() (i.e. all console output from growth() is hidden). Recruitment and summaries are done only once a year at the level of function fordyn().

Inspecting model outputs

Stand, species and cohort summaries

Among other outputs, function fordyn() calculates standard summary statistics that describe the structural and compositional state of the forest at each time step. For example, we can access stand-level statistics using:

fd$StandSummary
##   Step LeafAreaIndex TreeDensityLive TreeBasalAreaLive ShrubCoverLive MaxHeight
## 1    0      1.662671        552.0000          25.03330           3.75  800.0000
## 2    1      1.917676        546.3412          26.85181           3.75  824.6012
## 3    2      1.923576        540.7403          26.91888           3.75  828.9720
## 4    3      2.242156        535.1969          29.20546           3.75  855.8489

where we can observe an increase in stand basal area and leaf area index during years 1 and 3, but a reduction during the second (drier) year. Species-level analogous statistics are shown using:

fd$SpeciesSummary
##    Step Species LeafAreaIndex TreeDensityLive TreeBasalAreaLive ShrubCoverLive
## 1     0      54    0.81670117        168.0000         18.604545           0.00
## 2     0      65    0.04817502          0.0000          0.000000           3.75
## 3     0      68    0.79779523        384.0000          6.428754           0.00
## 4     1      54    0.91991321        166.2777         19.624708           0.00
## 5     1      65    0.04817502          0.0000          0.000000           3.75
## 6     1      68    0.94958748        380.0634          7.227105           0.00
## 7     2      54    0.91812052        164.5731         19.644870           0.00
## 8     2      65    0.04817502          0.0000          0.000000           3.75
## 9     2      68    0.95728001        376.1672          7.274006           0.00
## 10    3      54    1.02616274        162.8860         20.807347           0.00
## 11    3      65    0.04817502          0.0000          0.000000           3.75
## 12    3      68    1.16781834        372.3109          8.398109           0.00
##    MaxHeight
## 1   800.0000
## 2    80.0000
## 3   660.0000
## 4   824.6012
## 5    80.0000
## 6   685.5271
## 7   828.9720
## 8    80.0000
## 9   688.9028
## 10  855.8489
## 11   80.0000
## 12  722.4229

Tree/shrub tables

Another useful output of fordyn() are tables in long format with cohort structural information (i.e. DBH, height, density, etc) for each time step:

fd$TreeTable
##   Step Year Cohort Species             Name        N      DBH   Height Z50  Z95
## 1    0   NA  T1_54      54 Pinus halepensis 168.0000 37.55000 800.0000 750 3000
## 2    0   NA  T2_68      68     Quercus ilex 384.0000 14.60000 660.0000 750 3000
## 3    1 2001  T1_54      54 Pinus halepensis 166.2777 38.76498 824.6012 750 3000
## 4    1 2001  T2_68      68     Quercus ilex 380.0634 15.55999 685.5271 750 3000
## 5    2 2002  T1_54      54 Pinus halepensis 164.5731 38.98523 828.9720 750 3000
## 6    2 2002  T2_68      68     Quercus ilex 376.1672 15.69103 688.9028 750 3000
## 7    3 2003  T1_54      54 Pinus halepensis 162.8860 40.32938 855.8489 750 3000
## 8    3 2003  T2_68      68     Quercus ilex 372.3109 16.94700 722.4229 750 3000

The same can be shown for dead trees, where we see larger mortality rates during the second (drier) year:

fd$DeadTreeTable
##   Step Year Cohort Species             Name        N      DBH   Height Z50  Z95
## 1    1 2001  T1_54      54 Pinus halepensis 1.722253 38.76498 824.6012 750 3000
## 2    1 2001  T2_68      68     Quercus ilex 3.936579 15.55999 685.5271 750 3000
## 3    2 2002  T1_54      54 Pinus halepensis 1.704598 38.98523 828.9720 750 3000
## 4    2 2002  T2_68      68     Quercus ilex 3.896223 15.69103 688.9028 750 3000
## 5    3 2003  T1_54      54 Pinus halepensis 1.687123 40.32938 855.8489 750 3000
## 6    3 2003  T2_68      68     Quercus ilex 3.856281 16.94700 722.4229 750 3000

Accessing the output from function growth()

Since function fordyn() makes internal calls to function growth(), it stores the result in a vector called GrowthResults, which we can use to inspect intra-annual patterns of desired variables. For example, the following shows the leaf area for individuals of the three cohorts during the three consecutive years:

plot(fd$GrowthResults[[1]], "LeafArea", bySpecies = T)

plot(fd$GrowthResults[[2]], "LeafArea", bySpecies = T)

plot(fd$GrowthResults[[3]], "LeafArea", bySpecies = T)

where we see that the second (drier) year, resulted in a decrease in leaf area, due to the sink limitation caused by drought, and the third (normal) year lead to a recovery of leaf area. Note also the sink limitation to sapwood area (hence diameter) growth during the second year:

plot(fd$GrowthResults[[1]], "SapwoodArea", bySpecies = T)

plot(fd$GrowthResults[[2]], "SapwoodArea", bySpecies = T)

plot(fd$GrowthResults[[3]], "SapwoodArea", bySpecies = T)