Forest dynamics
Miquel De Caceres
2024-07-26
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():
examplesoil <- defaultSoilParams(4)
examplesoil## widths clay sand om nitrogen bd rfc
## 1 300 25 25 NA NA 1.5 25
## 2 700 25 25 NA NA 1.5 45
## 3 1000 25 25 NA NA 1.5 75
## 4 2000 25 25 NA NA 1.5 95As 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.3646 25.31498 812.5183 37.78378
## 3 550.7194 25.59839 824.9622 38.01864
## 4 550.0643 25.88293 837.3083 38.25412
## 5 549.3973 26.16814 849.5441 38.48999
## 6 548.7220 26.45403 861.6632 38.72609
## 7 548.0365 26.74035 873.6610 38.96234
## 8 547.3407 27.02699 885.5361 39.19868
## 9 546.6325 27.31373 897.2878 39.43508
## 10 545.9158 27.60076 908.9160 39.67150
## 11 545.1925 27.88817 920.4212 39.90795
## QuadraticMeanTreeDiameter HartBeckingIndex ShrubCoverLive BasalAreaDead
## 1 24.02949 53.20353 3.750000 0.00000000
## 2 24.17823 52.41401 3.112899 0.03949179
## 3 24.32743 51.65362 3.182481 0.04063946
## 4 24.47682 50.92229 3.253559 0.04181664
## 5 24.62625 50.21932 3.326337 0.04314122
## 6 24.77564 49.54346 3.400622 0.04425770
## 7 24.92493 48.89363 3.476457 0.04551891
## 8 25.07408 48.26862 3.553866 0.04680771
## 9 25.22306 47.66730 3.633780 0.04825718
## 10 25.37189 47.08834 3.713464 0.04947039
## 11 25.52056 46.53058 3.794796 0.05056136
## ShrubCoverDead BasalAreaCut ShrubCoverCut
## 1 0.000000000 0 0
## 2 0.005320281 0 0
## 3 0.004828445 0 0
## 4 0.004936335 0 0
## 5 0.005060626 0 0
## 6 0.005159621 0 0
## 7 0.005274778 0 0
## 8 0.005392333 0 0
## 9 0.005527954 0 0
## 10 0.005635651 0 0
## 11 0.005727148 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.6983 18.803102
## 5 1 Quercus coccifera 1 NA NA
## 6 1 Quercus ilex 1 383.6663 6.511881
## 7 2 Pinus halepensis 1 167.3916 19.002763
## 8 2 Quercus coccifera 1 NA NA
## 9 2 Quercus ilex 1 383.3277 6.595632
## 10 3 Pinus halepensis 1 167.0799 19.203066
## 11 3 Quercus coccifera 1 NA NA
## 12 3 Quercus ilex 1 382.9843 6.679863
## 13 4 Pinus halepensis 1 166.7622 19.403635
## 14 4 Quercus coccifera 1 NA NA
## 15 4 Quercus ilex 1 382.6351 6.764501
## 16 5 Pinus halepensis 1 166.4402 19.604490
## 17 5 Quercus coccifera 1 NA NA
## 18 5 Quercus ilex 1 382.2818 6.849543
## 19 6 Pinus halepensis 1 166.1130 19.805400
## 20 6 Quercus coccifera 1 NA NA
## 21 6 Quercus ilex 1 381.9235 6.934950
## 22 7 Pinus halepensis 1 165.7805 20.006279
## 23 7 Quercus coccifera 1 NA NA
## 24 7 Quercus ilex 1 381.5602 7.020708
## 25 8 Pinus halepensis 1 165.4418 20.206941
## 26 8 Quercus coccifera 1 NA NA
## 27 8 Quercus ilex 1 381.1907 7.106784
## 28 9 Pinus halepensis 1 165.0987 20.407551
## 29 9 Quercus coccifera 1 NA NA
## 30 9 Quercus ilex 1 380.8171 7.193207
## 31 10 Pinus halepensis 1 164.7521 20.608184
## 32 10 Quercus coccifera 1 NA NA
## 33 10 Quercus ilex 1 380.4404 7.279990
## 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.033827461 NA 0 NA
## 5 3.112899 NA 0.005320281 NA 0
## 6 NA 0.005664332 NA 0 NA
## 7 NA 0.034814617 NA 0 NA
## 8 3.182481 NA 0.004828445 NA 0
## 9 NA 0.005824846 NA 0 NA
## 10 NA 0.035827331 NA 0 NA
## 11 3.253559 NA 0.004936335 NA 0
## 12 NA 0.005989313 NA 0 NA
## 13 NA 0.036966589 NA 0 NA
## 14 3.326337 NA 0.005060626 NA 0
## 15 NA 0.006174631 NA 0 NA
## 16 NA 0.037927725 NA 0 NA
## 17 3.400622 NA 0.005159621 NA 0
## 18 NA 0.006329972 NA 0 NA
## 19 NA 0.039013042 NA 0 NA
## 20 3.476457 NA 0.005274778 NA 0
## 21 NA 0.006505865 NA 0 NA
## 22 NA 0.040122162 NA 0 NA
## 23 3.553866 NA 0.005392333 NA 0
## 24 NA 0.006685549 NA 0 NA
## 25 NA 0.041369169 NA 0 NA
## 26 3.633780 NA 0.005527954 NA 0
## 27 NA 0.006888013 NA 0 NA
## 28 NA 0.042413749 NA 0 NA
## 29 3.713464 NA 0.005635651 NA 0
## 30 NA 0.007056640 NA 0 NA
## 31 NA 0.043353544 NA 0 NA
## 32 3.794796 NA 0.005727148 NA 0
## 33 NA 0.007207812 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.78378 812.5183 167.6983 100 600
## 4 1 2001 T2_168 Quercus ilex 14.70048 663.1356 383.6663 300 1000
## 5 2 2002 T1_148 Pinus halepensis 38.01864 824.9622 167.3916 100 600
## 6 2 2002 T2_168 Quercus ilex 14.80124 666.2859 383.3277 300 1000
## 7 3 2003 T1_148 Pinus halepensis 38.25412 837.3083 167.0799 100 600
## 8 3 2003 T2_168 Quercus ilex 14.90213 669.4451 382.9843 300 1000
## 9 4 2004 T1_148 Pinus halepensis 38.48999 849.5441 166.7622 100 600
## 10 4 2004 T2_168 Quercus ilex 15.00308 672.6104 382.6351 300 1000
## 11 5 2005 T1_148 Pinus halepensis 38.72609 861.6632 166.4402 100 600
## 12 5 2005 T2_168 Quercus ilex 15.10407 675.7801 382.2818 300 1000
## 13 6 2006 T1_148 Pinus halepensis 38.96234 873.6610 166.1130 100 600
## 14 6 2006 T2_168 Quercus ilex 15.20507 678.9523 381.9235 300 1000
## 15 7 2007 T1_148 Pinus halepensis 39.19868 885.5361 165.7805 100 600
## 16 7 2007 T2_168 Quercus ilex 15.30608 682.1260 381.5602 300 1000
## 17 8 2008 T1_148 Pinus halepensis 39.43508 897.2878 165.4418 100 600
## 18 8 2008 T2_168 Quercus ilex 15.40708 685.3002 381.1907 300 1000
## 19 9 2009 T1_148 Pinus halepensis 39.67150 908.9160 165.0987 100 600
## 20 9 2009 T2_168 Quercus ilex 15.50808 688.4738 380.8171 300 1000
## 21 10 2010 T1_148 Pinus halepensis 39.90795 920.4212 164.7521 100 600
## 22 10 2010 T2_168 Quercus ilex 15.60907 691.6461 380.4404 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.78378 812.5183 0.3016953 0
## 2 1 2001 T2_168 Quercus ilex 14.70048 663.1356 0.3337304 0
## 3 2 2002 T1_148 Pinus halepensis 38.01864 824.9622 0.3066752 0
## 4 2 2002 T2_168 Quercus ilex 14.80124 666.2859 0.3385309 0
## 5 3 2003 T1_148 Pinus halepensis 38.25412 837.3083 0.3117225 0
## 6 3 2003 T2_168 Quercus ilex 14.90213 669.4451 0.3433922 0
## 7 4 2004 T1_148 Pinus halepensis 38.48999 849.5441 0.3177049 0
## 8 4 2004 T2_168 Quercus ilex 15.00308 672.6104 0.3492690 0
## 9 5 2005 T1_148 Pinus halepensis 38.72609 861.6632 0.3220027 0
## 10 5 2005 T2_168 Quercus ilex 15.10407 675.7801 0.3532839 0
## 11 6 2006 T1_148 Pinus halepensis 38.96234 873.6610 0.3272124 0
## 12 6 2006 T2_168 Quercus ilex 15.20507 678.9523 0.3582928 0
## 13 7 2007 T1_148 Pinus halepensis 39.19868 885.5361 0.3324693 0
## 14 7 2007 T2_168 Quercus ilex 15.30608 682.1260 0.3633450 0
## 15 8 2008 T1_148 Pinus halepensis 39.43508 897.2878 0.3387049 0
## 16 8 2008 T2_168 Quercus ilex 15.40708 685.3002 0.3694563 0
## 17 9 2009 T1_148 Pinus halepensis 39.67150 908.9160 0.3431305 0
## 18 9 2009 T2_168 Quercus ilex 15.50808 688.4738 0.3735871 0
## 19 10 2010 T1_148 Pinus halepensis 39.90795 920.4212 0.3465898 0
## 20 10 2010 T2_168 Quercus ilex 15.60907 691.6461 0.3766685 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.78378 812.5183 9.655462 100 600
## 2 1 2001 T2_168 Quercus ilex 14.70048 663.1356 383.666270 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
## 1 0 NA T1_148 Pinus halepensis 37.550000 800.00000 168.0000 100.0000
## 2 0 NA T2_168 Quercus ilex 14.600000 660.00000 384.0000 300.0000
## 3 1 2001 T1_148 Pinus halepensis 37.783784 812.51834 158.0428 100.0000
## 4 1 2001 T3_168 Quercus ilex 1.000000 47.23629 3000.0000 300.0000
## 5 2 2002 T1_148 Pinus halepensis 38.019501 824.89925 157.8702 100.0000
## 6 2 2002 T3_168 Quercus ilex 1.101965 53.48699 2705.7850 300.0000
## 7 3 2003 T1_148 Pinus halepensis 38.255747 837.16631 157.6947 100.0000
## 8 3 2003 T3_168 Quercus ilex 1.205074 59.77879 2460.3536 300.0000
## 9 4 2004 T1_148 Pinus halepensis 38.492386 849.32474 157.5159 100.0000
## 10 4 2004 T3_168 Quercus ilex 1.308004 66.06134 2255.0521 300.0000
## 11 5 2005 T1_148 Pinus halepensis 38.729290 861.36881 157.3347 100.0000
## 12 5 2005 T2_168 Quercus ilex 1.350610 70.03350 2501.4760 281.5475
## 13 6 2006 T1_148 Pinus halepensis 38.963966 873.17410 157.1508 100.0000
## 14 6 2006 T2_168 Quercus ilex 1.453169 76.29847 2016.3537 281.5475
## 15 7 2007 T1_148 Pinus halepensis 39.199600 884.90308 156.9641 100.0000
## 16 7 2007 T2_168 Quercus ilex 1.555897 82.57401 1875.1291 281.5475
## 17 8 2008 T1_148 Pinus halepensis 39.435830 896.53785 156.7738 100.0000
## 18 8 2008 T2_168 Quercus ilex 1.658767 88.86540 1751.7500 281.5475
## 19 9 2009 T1_148 Pinus halepensis 39.672446 908.06840 156.5811 100.0000
## 20 9 2009 T2_168 Quercus ilex 1.761721 95.16919 1643.1273 281.5475
## 21 10 2010 T1_148 Pinus halepensis 39.909316 919.48909 156.3865 100.0000
## 22 10 2010 T2_168 Quercus ilex 1.864722 101.48309 1546.8185 281.5475
## Z95
## 1 600
## 2 1000
## 3 600
## 4 1000
## 5 600
## 6 1000
## 7 600
## 8 1000
## 9 600
## 10 1000
## 11 600
## 12 1000
## 13 600
## 14 1000
## 15 600
## 16 1000
## 17 600
## 18 1000
## 19 600
## 20 1000
## 21 600
## 22 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).