Forest dynamics
Miquel De Caceres
2024-04-16
Source:vignettes/runmodels/ForestDynamics.Rmd
ForestDynamics.RmdAbout this vignette
This document describes how to run the forest dynamics model of
medfate, described in De Cáceres et al. (2023) and
implemented in function fordyn(). This document is meant to
teach users to run the simulation model with function
fordyn(). Details of the model design and formulation can
be found at the corresponding chapters of the medfate
book.
Because the model builds on the growth and water balance models, the
reader is assumed here to be familiarized with spwb() and
growth() (otherwise read vignettes Basic
water balance and Forest
growth).
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().
Model inputs are explained in greater detail in vignettes Understanding
model inputs and Preparing
model 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(4)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(exampleforest)
exampleforest## $treeData
## Species N DBH Height Z50 Z95
## 1 Pinus halepensis 168 37.55 800 100 600
## 2 Quercus ilex 384 14.60 660 300 1000
##
## $shrubData
## Species Cover Height Z50 Z95
## 1 Quercus coccifera 3.75 80 200 1000
##
## $herbCover
## [1] 10
##
## $herbHeight
## [1] 20
##
## $seedBank
## [1] Species Percent
## <0 rows> (or 0-length row.names)
##
## attr(,"class")
## [1] "forest" "list"Importantly, a data frame with daily weather for the period to be simulated is required. Here we use the default data frame included with the package:
## dates MinTemperature MaxTemperature Precipitation MinRelativeHumidity
## 1 2001-01-01 -0.5934215 6.287950 4.869109 65.15411
## 2 2001-01-02 -2.3662458 4.569737 2.498292 57.43761
## 3 2001-01-03 -3.8541036 2.661951 0.000000 58.77432
## 4 2001-01-04 -1.8744860 3.097705 5.796973 66.84256
## 5 2001-01-05 0.3288287 7.551532 1.884401 62.97656
## 6 2001-01-06 0.5461322 7.186784 13.359801 74.25754
## MaxRelativeHumidity Radiation WindSpeed
## 1 100.00000 12.89251 2.000000
## 2 94.71780 13.03079 7.662544
## 3 94.66823 16.90722 2.000000
## 4 95.80950 11.07275 2.000000
## 5 100.00000 13.45205 7.581347
## 6 100.00000 12.84841 6.570501Finally, 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 ten-year weather input by repeating the example weather data frame ten times.
meteo <- rbind(examplemeteo, examplemeteo, examplemeteo, examplemeteo,
examplemeteo, examplemeteo, examplemeteo, examplemeteo,
examplemeteo, examplemeteo)
meteo$dates = seq(as.Date("2001-01-01"),
as.Date("2010-12-29"), by="day")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(exampleforest, examplesoil, SpParamsMED, meteo, control,
latitude = 41.82592, elevation = 100)## Simulating year 2001 (1/10): (a) Growth/mortality## Package 'meteoland' [ver. 2.2.1]## , (b) Regeneration nT = 2 nS = 1
## Simulating year 2002 (2/10): (a) Growth/mortality, (b) Regeneration nT = 2 nS = 1
## Simulating year 2003 (3/10): (a) Growth/mortality, (b) Regeneration nT = 2 nS = 1
## Simulating year 2004 (4/10): (a) Growth/mortality, (b) Regeneration nT = 2 nS = 1
## Simulating year 2005 (5/10): (a) Growth/mortality, (b) Regeneration nT = 2 nS = 1
## Simulating year 2006 (6/10): (a) Growth/mortality, (b) Regeneration nT = 2 nS = 1
## Simulating year 2007 (7/10): (a) Growth/mortality, (b) Regeneration nT = 2 nS = 1
## Simulating year 2008 (8/10): (a) Growth/mortality, (b) Regeneration nT = 2 nS = 1
## Simulating year 2009 (9/10): (a) Growth/mortality, (b) Regeneration nT = 2 nS = 1
## Simulating year 2010 (10/10): (a) Growth/mortality, (b) Regeneration nT = 2 nS = 1It 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 and plots
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 NumTreeSpecies NumTreeCohorts NumShrubSpecies NumShrubCohorts
## 1 0 2 2 1 1
## 2 1 2 2 1 1
## 3 2 2 2 1 1
## 4 3 2 2 1 1
## 5 4 2 2 1 1
## 6 5 2 2 1 1
## 7 6 2 2 1 1
## 8 7 2 2 1 1
## 9 8 2 2 1 1
## 10 9 2 2 1 1
## 11 10 2 2 1 1
## TreeDensityLive TreeBasalAreaLive DominantTreeHeight DominantTreeDiameter
## 1 552.0000 25.03330 800.0000 37.55000
## 2 551.3644 25.32558 812.9229 37.79138
## 3 550.7186 25.61949 825.7531 38.03365
## 4 550.0625 25.91455 838.4733 38.27647
## 5 549.3941 26.21037 851.0746 38.51966
## 6 548.7170 26.50695 863.5531 38.76313
## 7 548.0292 26.80401 875.9057 39.00681
## 8 547.3306 27.10150 888.1313 39.25065
## 9 546.6193 27.39923 900.2293 39.49462
## 10 545.8989 27.69736 912.1994 39.73869
## 11 545.1715 27.99601 924.0417 39.98284
## QuadraticMeanTreeDiameter HartBeckingIndex ShrubCoverLive BasalAreaDead
## 1 24.02949 53.20353 3.750000 0.00000000
## 2 24.18329 52.38794 3.114844 0.03952091
## 3 24.33747 51.60419 3.186505 0.04071233
## 4 24.49181 50.85161 3.259778 0.04193504
## 5 24.64618 50.12915 3.334870 0.04330784
## 6 24.80052 49.43525 3.411581 0.04447395
## 7 24.95474 48.76866 3.489955 0.04578760
## 8 25.10885 48.12801 3.570021 0.04713139
## 9 25.26282 47.51212 3.652743 0.04863972
## 10 25.41664 46.91958 3.735300 0.04991253
## 11 25.57034 46.34916 3.819628 0.05106411
## ShrubCoverDead BasalAreaCut ShrubCoverCut
## 1 0.000000000 0 0
## 2 0.005321369 0 0
## 3 0.004832600 0 0
## 4 0.004943751 0 0
## 5 0.005071520 0 0
## 6 0.005174105 0 0
## 7 0.005293056 0 0
## 8 0.005414582 0 0
## 9 0.005554446 0 0
## 10 0.005666425 0 0
## 11 0.005762201 0 0Species-level analogous statistics are shown using:
fd$SpeciesSummary## Step Species NumCohorts TreeDensityLive TreeBasalAreaLive
## 1 0 Pinus halepensis 1 168.0000 18.604547
## 2 0 Quercus coccifera 1 NA NA
## 3 0 Quercus ilex 1 384.0000 6.428755
## 4 1 Pinus halepensis 1 167.6982 18.810652
## 5 1 Quercus coccifera 1 NA NA
## 6 1 Quercus ilex 1 383.6662 6.514925
## 7 2 Pinus halepensis 1 167.3912 19.017722
## 8 2 Quercus coccifera 1 NA NA
## 9 2 Quercus ilex 1 383.3274 6.601763
## 10 3 Pinus halepensis 1 167.0790 19.225402
## 11 3 Quercus coccifera 1 NA NA
## 12 3 Quercus ilex 1 382.9835 6.689144
## 13 4 Pinus halepensis 1 166.7605 19.433371
## 14 4 Quercus coccifera 1 NA NA
## 15 4 Quercus ilex 1 382.6335 6.776996
## 16 5 Pinus halepensis 1 166.4376 19.641699
## 17 5 Quercus coccifera 1 NA NA
## 18 5 Quercus ilex 1 382.2794 6.865246
## 19 6 Pinus halepensis 1 166.1092 19.850176
## 20 6 Quercus coccifera 1 NA NA
## 21 6 Quercus ilex 1 381.9200 6.953837
## 22 7 Pinus halepensis 1 165.7753 20.058727
## 23 7 Quercus coccifera 1 NA NA
## 24 7 Quercus ilex 1 381.5554 7.042768
## 25 8 Pinus halepensis 1 165.4349 20.267161
## 26 8 Quercus coccifera 1 NA NA
## 27 8 Quercus ilex 1 381.1844 7.132067
## 28 9 Pinus halepensis 1 165.0899 20.475637
## 29 9 Quercus coccifera 1 NA NA
## 30 9 Quercus ilex 1 380.8090 7.221723
## 31 10 Pinus halepensis 1 164.7411 20.684227
## 32 10 Quercus coccifera 1 NA NA
## 33 10 Quercus ilex 1 380.4304 7.311778
## ShrubCoverLive BasalAreaDead ShrubCoverDead BasalAreaCut ShrubCoverCut
## 1 NA 0.000000000 NA 0 NA
## 2 3.750000 NA 0.000000000 NA 0
## 3 NA 0.000000000 NA 0 NA
## 4 NA 0.033852351 NA 0 NA
## 5 3.114844 NA 0.005321369 NA 0
## 6 NA 0.005668555 NA 0 NA
## 7 NA 0.034877172 NA 0 NA
## 8 3.186505 NA 0.004832600 NA 0
## 9 NA 0.005835159 NA 0 NA
## 10 NA 0.035929045 NA 0 NA
## 11 3.259778 NA 0.004943751 NA 0
## 12 NA 0.006005997 NA 0 NA
## 13 NA 0.037109757 NA 0 NA
## 14 3.334870 NA 0.005071520 NA 0
## 15 NA 0.006198079 NA 0 NA
## 16 NA 0.038113617 NA 0 NA
## 17 3.411581 NA 0.005174105 NA 0
## 18 NA 0.006360328 NA 0 NA
## 19 NA 0.039244150 NA 0 NA
## 20 3.489955 NA 0.005293056 NA 0
## 21 NA 0.006543447 NA 0 NA
## 22 NA 0.040400743 NA 0 NA
## 23 3.570021 NA 0.005414582 NA 0
## 24 NA 0.006730642 NA 0 NA
## 25 NA 0.041698592 NA 0 NA
## 26 3.652743 NA 0.005554446 NA 0
## 27 NA 0.006941133 NA 0 NA
## 28 NA 0.042794702 NA 0 NA
## 29 3.735300 NA 0.005666425 NA 0
## 30 NA 0.007117832 NA 0 NA
## 31 NA 0.043786924 NA 0 NA
## 32 3.819628 NA 0.005762201 NA 0
## 33 NA 0.007277189 NA 0 NAPackage medfate provides a simple plot
function for objects of class fordyn. For example, we can
show the interannual variation in stand-level basal area using:
plot(fd, type = "StandBasalArea")
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 DBH Height N Z50 Z95
## 1 0 NA T1_148 Pinus halepensis 37.55000 800.0000 168.0000 100 600
## 2 0 NA T2_168 Quercus ilex 14.60000 660.0000 384.0000 300 1000
## 3 1 2001 T1_148 Pinus halepensis 37.79138 812.9229 167.6982 100 600
## 4 1 2001 T2_168 Quercus ilex 14.70392 663.2432 383.6662 300 1000
## 5 2 2002 T1_148 Pinus halepensis 38.03365 825.7531 167.3912 100 600
## 6 2 2002 T2_168 Quercus ilex 14.80813 666.5017 383.3274 300 1000
## 7 3 2003 T1_148 Pinus halepensis 38.27647 838.4733 167.0790 100 600
## 8 3 2003 T2_168 Quercus ilex 14.91249 669.7703 382.9835 300 1000
## 9 4 2004 T1_148 Pinus halepensis 38.51966 851.0746 166.7605 100 600
## 10 4 2004 T2_168 Quercus ilex 15.01696 673.0463 382.6335 300 1000
## 11 5 2005 T1_148 Pinus halepensis 38.76313 863.5531 166.4376 100 600
## 12 5 2005 T2_168 Quercus ilex 15.12142 676.3253 382.2794 300 1000
## 13 6 2006 T1_148 Pinus halepensis 39.00681 875.9057 166.1092 100 600
## 14 6 2006 T2_168 Quercus ilex 15.22583 679.6049 381.9200 300 1000
## 15 7 2007 T1_148 Pinus halepensis 39.25065 888.1313 165.7753 100 600
## 16 7 2007 T2_168 Quercus ilex 15.33020 682.8846 381.5554 300 1000
## 17 8 2008 T1_148 Pinus halepensis 39.49462 900.2293 165.4349 100 600
## 18 8 2008 T2_168 Quercus ilex 15.43459 686.1652 381.1844 300 1000
## 19 9 2009 T1_148 Pinus halepensis 39.73869 912.1994 165.0899 100 600
## 20 9 2009 T2_168 Quercus ilex 15.53896 689.4444 380.8090 300 1000
## 21 10 2010 T1_148 Pinus halepensis 39.98284 924.0417 164.7411 100 600
## 22 10 2010 T2_168 Quercus ilex 15.64332 692.7223 380.4304 300 1000The same can be shown for dead trees:
fd$DeadTreeTable## Step Year Cohort Species DBH Height N N_starvation
## 1 1 2001 T1_148 Pinus halepensis 37.79138 812.9229 0.3017959 0
## 2 1 2001 T2_168 Quercus ilex 14.70392 663.2432 0.3338232 0
## 3 2 2002 T1_148 Pinus halepensis 38.03365 825.7531 0.3069838 0
## 4 2 2002 T2_168 Quercus ilex 14.80813 666.5017 0.3388150 0
## 5 3 2003 T1_148 Pinus halepensis 38.27647 838.4733 0.3122425 0
## 6 3 2003 T2_168 Quercus ilex 14.91249 669.7703 0.3438703 0
## 7 4 2004 T1_148 Pinus halepensis 38.51966 851.0746 0.3184441 0
## 8 4 2004 T2_168 Quercus ilex 15.01696 673.0463 0.3499475 0
## 9 5 2005 T1_148 Pinus halepensis 38.76313 863.5531 0.3229628 0
## 10 5 2005 T2_168 Quercus ilex 15.12142 676.3253 0.3541639 0
## 11 6 2006 T1_148 Pinus halepensis 39.00681 875.9057 0.3284008 0
## 12 6 2006 T2_168 Quercus ilex 15.22583 679.6049 0.3593805 0
## 13 7 2007 T1_148 Pinus halepensis 39.25065 888.1313 0.3338918 0
## 14 7 2007 T2_168 Quercus ilex 15.33020 682.8846 0.3646453 0
## 15 8 2008 T1_148 Pinus halepensis 39.49462 900.2293 0.3403734 0
## 16 8 2008 T2_168 Quercus ilex 15.43459 686.1652 0.3709796 0
## 17 9 2009 T1_148 Pinus halepensis 39.73869 912.1994 0.3450428 0
## 18 9 2009 T2_168 Quercus ilex 15.53896 689.4444 0.3753307 0
## 19 10 2010 T1_148 Pinus halepensis 39.98284 924.0417 0.3487443 0
## 20 10 2010 T2_168 Quercus ilex 15.64332 692.7223 0.3786308 0
## N_dessication N_burnt Z50 Z95
## 1 0 0 100 600
## 2 0 0 300 1000
## 3 0 0 100 600
## 4 0 0 300 1000
## 5 0 0 100 600
## 6 0 0 300 1000
## 7 0 0 100 600
## 8 0 0 300 1000
## 9 0 0 100 600
## 10 0 0 300 1000
## 11 0 0 100 600
## 12 0 0 300 1000
## 13 0 0 100 600
## 14 0 0 300 1000
## 15 0 0 100 600
## 16 0 0 300 1000
## 17 0 0 100 600
## 18 0 0 300 1000
## 19 0 0 100 600
## 20 0 0 300 1000Accessing 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 second year:
plot(fd$GrowthResults[[2]], "LeafArea", bySpecies = T)
Instead of examining year by year, it is possible to plot the whole
series of results by passing a fordyn object to the
plot() function:
plot(fd, "LeafArea")
Finally, we can create interactive plots for particular steps using
function shinyplot(), e.g.:
shinyplot(fd$GrowthResults[[1]])Forest dynamics including management
The package allows including forest management in simulations of
forest dynamics. This is done in a very flexible manner, in the sense
that fordyn() allows the user to supply an arbitrary
function implementing a desired management strategy for the stand whose
dynamics are to be simulated. The package includes, however, an in-built
default function called defaultManagementFunction() along
with a flexible parameterization, a list with defaults provided by
function defaultManagementArguments().
Here we provide an example of simulations including forest management:
# Default arguments
args <- defaultManagementArguments()
# Here one can modify defaults before calling fordyn()
#
# Simulation
fd<-fordyn(exampleforest, examplesoil, SpParamsMED, meteo, control,
latitude = 41.82592, elevation = 100,
management_function = defaultManagementFunction,
management_args = args)## Simulating year 2001 (1/10): (a) Growth/mortality & management [thinning], (b) Regeneration nT = 2 nS = 2
## Simulating year 2002 (2/10): (a) Growth/mortality & management [none], (b) Regeneration nT = 2 nS = 2
## Simulating year 2003 (3/10): (a) Growth/mortality & management [none], (b) Regeneration nT = 2 nS = 2
## Simulating year 2004 (4/10): (a) Growth/mortality & management [none], (b) Regeneration nT = 2 nS = 2
## Simulating year 2005 (5/10): (a) Growth/mortality & management [none], (b) Regeneration nT = 2 nS = 2
## Simulating year 2006 (6/10): (a) Growth/mortality & management [none], (b) Regeneration nT = 2 nS = 2
## Simulating year 2007 (7/10): (a) Growth/mortality & management [none], (b) Regeneration nT = 2 nS = 2
## Simulating year 2008 (8/10): (a) Growth/mortality & management [none], (b) Regeneration nT = 2 nS = 2
## Simulating year 2009 (9/10): (a) Growth/mortality & management [none], (b) Regeneration nT = 2 nS = 2
## Simulating year 2010 (10/10): (a) Growth/mortality & management [none], (b) Regeneration nT = 2 nS = 2When management is included in simulations, two additional tables are produced, corresponding to the trees and shrubs that were cut, e.g.:
fd$CutTreeTable## Step Year Cohort Species DBH Height N Z50 Z95
## 1 1 2001 T1_148 Pinus halepensis 37.79138 812.9229 9.652777 100 600
## 2 1 2001 T2_168 Quercus ilex 14.70392 663.2432 383.666177 300 1000Management parameters were those of an irregular model with thinning interventions from ‘below’, indicating that smaller trees were to be cut earlier:
args$type## [1] "irregular"
args$thinning## [1] "below"Note that in this example, there is resprouting of Quercus ilex after the thinning intervention, evidenced by the new cohort (T3_168) appearing in year 2001:
fd$TreeTable## Step Year Cohort Species DBH Height N Z50 Z95
## 1 0 NA T1_148 Pinus halepensis 37.550000 800.00000 168.0000 100 600
## 2 0 NA T2_168 Quercus ilex 14.600000 660.00000 384.0000 300 1000
## 3 1 2001 T1_148 Pinus halepensis 37.791380 812.92290 158.0454 100 600
## 4 1 2001 T3_168 Quercus ilex 1.000000 47.23629 3000.0000 300 1000
## 5 2 2002 T1_148 Pinus halepensis 38.034209 825.67126 157.8726 100 600
## 6 2 2002 T3_168 Quercus ilex 1.105169 53.68427 2697.4474 300 1000
## 7 3 2003 T1_148 Pinus halepensis 38.277486 838.29312 157.6969 100 600
## 8 3 2003 T3_168 Quercus ilex 1.210647 60.12140 2448.3140 300 1000
## 9 4 2004 T1_148 Pinus halepensis 38.521106 850.79615 157.5177 100 600
## 10 4 2004 T3_168 Quercus ilex 1.315957 66.54987 2240.5659 300 1000
## 11 5 2005 T1_148 Pinus halepensis 38.764979 863.17667 157.3361 100 600
## 12 5 2005 T3_168 Quercus ilex 1.421362 72.98707 2064.3551 300 1000
## 13 6 2006 T1_148 Pinus halepensis 39.009050 875.43259 157.1515 100 600
## 14 6 2006 T3_168 Quercus ilex 1.526941 79.43806 1912.9535 300 1000
## 15 7 2007 T1_148 Pinus halepensis 39.253278 887.56281 156.9638 100 600
## 16 7 2007 T3_168 Quercus ilex 1.632625 85.89959 1781.5867 300 1000
## 17 8 2008 T1_148 Pinus halepensis 39.497634 899.56701 156.7726 100 600
## 18 8 2008 T3_168 Quercus ilex 1.738347 92.36876 1666.6236 300 1000
## 19 9 2009 T1_148 Pinus halepensis 39.742101 911.44543 156.5787 100 600
## 20 9 2009 T3_168 Quercus ilex 1.844205 98.85374 1565.1133 300 1000
## 21 10 2010 T1_148 Pinus halepensis 39.986660 923.19825 156.3829 100 600
## 22 10 2010 T3_168 Quercus ilex 1.950250 105.35771 1474.8061 300 1000References
- De Cáceres M, Molowny-Horas R, Cabon A, Martínez-Vilalta J, Mencuccini M, García-Valdés R, Nadal-Sala D, Sabaté S, Martin-StPaul N, Morin X, D’Adamo F, Batllori E, Améztegui A (2023) MEDFATE 2.9.3: A trait-enabled model to simulate Mediterranean forest function and dynamics at regional scales. Geoscientific Model Development 16: 3165-3201 (https://doi.org/10.5194/gmd-16-3165-2023).