Model description, parameterization and evaluation of formind for tropical montane cloud forest in central Veracruz, Mexico

1. 
   FORMIND’s submodels for recruitment, mortality, light competition, growth, and logging
   
2. 
   the model parameterization for tropical montane cloud forest in central Veracruz, Mexico,
   
3. 
   the model evaluation,
   
4. 
   and a list of tree species in the study area (Table S3).
   

1. DESCRIPTION OF SUBMODELS OF FORMIND

If irradiance at the forest floor (I_floor) matches light requirements for establishment of trees of a given PFT (I_min, I_max), that is,

(1) Basic mortality: Each tree has a basic PFT-specific probability to die (m_B).

(2) Size-dependent mortality: Small trees experience additional mortality (m_D) depending on their actual dbh (D),

where m_max is the maximum size-dependent mortality of small trees and D_mort is the dbh up to which mortality is increased. Basic and size-dependent mortality are added,

For cohorts with more than 100 individuals and a dbh less than 10 cm, the number of dying trees (N_d) is calculated as

where N_t is the number of trees in the cohort; for the non-integral part of N_d it is stochastically determined if an additional tree dies. Likewise, for cohorts with less than 100 individuals or dbh of 10 cm or greater it is stochastically determined if a tree dies.

(3) Self-thinning: In dense patches mortality is increased due to competition for space among the trees. When the sum of the crown area of all trees (which have their crown in this height layer) in a given height layer exceeds the patch area, trees are randomly removed until the crowns of all trees “fit” into the patch.

(4) Gap creation: Large falling trees kill a proportion of the trees in the patch where their crown hits the ground. Trees with a dbh greater than D_fall fall over with the probability p_fall. The probability that a tree in the target patch is killed (p_k) is proportional to the ratio between the crown projection area of the falling tree (CA) and patch size (a),

Again, for cohorts with more than 100 individuals and a dbh less than 10 cm, the number of dying trees is calculated deterministically, otherwise it is stochastically determined if a tree dies. Only trees that do not overtop the falling tree by more than 1 m can be killed.

The vertical distribution of leaf area determines the light climate in each patch. This is accounted for by dividing each patch into vertical layers of the width Δh. The crowns of most trees span over several height layers, and thus contribute to the leaf area index (LAI) in each height layer i that contains a part of the tree crown (),

with N_t being the number of trees, CA_t crown area, CL_t crown depth, LAI_t LAI of cohort t, and a patch area. When crowns of large trees exceed the patch area, “overhanging” crown and leaf area is distributed evenly to the corresponding height layers of the four neighboring patches.

The cumulative LAI () adds up the LAI in all height layers above layer i,

.

The available light at the crown of each tree (I_t) is calculated via an extinction law

with I₀ being the average irradiance above the forest canopy, LAIC_t the cumulative LAI in all height layers above the crown of the tree and k the light extinction coefficient of the forest.

The calculation of light extinction within the forest canopy and leaf-level rates of photosynthesis follows Thornley and Johnson (1990) and Monsi and Saeki (1953). Both, incident irradiance and rate of photosynthesis need to be considered on the level of single leaves (that is, per unit leaf area) and the level of an entire tree crown (per unit crown projection area), denoted by l and t respectively. The rate of single-leaf photosynthesis (P_l) is modelled as a saturating function of the incident light on the leaf (I_l),

with p_max being the maximum rate of photosynthesis and α the initial slope of the light-response curve. The irradiance incident on the surface of a leaf within the canopy of a tree is

with I_t being the irradiance incident on the tree crown, k the light extinction coefficient of the forest, and m the light transmission coefficient of leaves. Self-shading of the tree canopy is accounted for by integrating P_l over the LAI of tree t (LAI_t) and the resulting instantaneous rate of photosynthesis of the tree (P_t, per unit crown projection area) is

where L is the cumulative LAI of the tree. Solving the integral leads to

(Thornley and Johnson 1990).

To calculate total annual gross biomass production of the tree (PB), the rate of photosynthesis is multiplied by the length of the photosynthetic active period per year (S), the crown area of the tree (CA) and a conversion coefficient (codm) from absorbed CO₂ to organic dry mass,

.

S is calculated as

,

where s_d is the average daily photosynthetic active period.

Respiration processes can be divided into growth respiration during the build-up of new biomass (parameter r_g) and maintenance respiration of living biomass (parameters r₀, r₁, r₂). Growth respiration is assumed to be a fixed fraction of net biomass production, whereas maintenance respiration is assumed to be non-linearly dependent on the living biomass (B) of the tree. Thus, biomass increment B_inc of the tree is calculated as

.

Maintenance respiration parameters are fitted such that measured maximum diameter increment values for each PFT are reproduced. Additionally it is assured that. The new biomass () is then translated into the new dbh of the tree via a table function. The table is filled once at the beginning of a simulation for each PFT and assigns each possible dbh (in steps of one mm) the corresponding biomass on the basis of the equation

with D being tree dbh, H being the height of the tree, f the form factor that corrects the deviation of the stem from the idealised conical shape, ρ the wood density and sw the fraction of stem wood biomass from total tree biomass. Between tabulated values linear interpolation is applied.

From the new dbh (sometimes converted into cm) all other variables describing the geometry of the tree are derived. Height (H) is calculated as

.

Crown depth (CL) is a constant fraction of height,

Crown diameter (CD) is assumed to be proportional to the stem diameter,

The crown is assumed to be a cylinder, hence crown area (CA) is

LAI of a tree (LAI) is fixed to L. Stem volume (SV) is

Bole volume (stem volume below the crown, BV) is derived from geometrical properties of the frustum of a cone (Bronstein and Semendjajew 1991; Ditzer 1999)

with x being

Selective logging: The program keeps track of harvestable trees that comply with defined criteria of the logging scenarios (for example, commercial PFTs, minimum and maximum allowed dbh thresholds for harvesting). Before the logging module is applied, it is evaluated whether the minimum criterion (that is, minimum number of trees to be extracted per hectare) can be fulfilled taking potential logging damage into account. If the minimum criterion is met, a logging operation takes place; otherwise logging is omitted in the respective hectare. Patches are visited randomly and the largest harvestable tree of the patch is logged, until all patches have been visited at least once or the harvest target has been met. Then, patches are visited randomly until the harvest target is met. Reduced impact logging is assumed. This means, that falling logged trees are directed to previously damaged patches, if possible, to reduce logging damages by falling trees. Vegetation damage in the patch where the crown of the logged tree hits the ground is simulated in the same way as described for naturally falling trees (see Mortality (4) Gap creation), except that harvestable trees in the target patch are prevented from being damaged. When a tree is directed to a previously damaged patch, only the difference between potential damage caused by its crown and the previously damaged proportion of the trees is damaged. Additional logging damage due to skidding apply to the entire hectare, independently of the location of logged trees. The intensity of this global damage is defined for each logging scenario depending on the logging intensity.

2. MODEL PARAMETERS
Model parameters for tropical montane cloud forest (TMCF) in central Veracruz, Mexico, were obtained for environment, tree geometry, tree growth, recruitment, mortality, and technical issues (Table S2). In the description of model parameters we focus on biologically relevant parameters, more technical details can be found in Köhler (2000).
Environmental parameters

The climate of central Veracruz is mild and humid and we assume that trees can grow throughout the whole year. Mean annual light intensity above the canopy was estimated to be 600 µmol(photons)·m^-2s^-1. Hafkenscheid (2000) measured average values of 1150 µmol(photons)·m^-2s^-1 on a sunny day and 260 µmol(photons)·m^-2s^-1 on a totally overcast day in the Jamaican Blue Mountains. As annual average light intensity he calculated 650 µmol(photons)·m^-2s^-1. Measurements of the light extinction coefficient for TMCF are not available. In tall tropical lowland forests a commonly observed value is 0.7 (Kira 1978). Hafkenscheid (2000) estimated the light extinction coefficient for short-statured upper montane rain forest to be 0.5. Hence, a light extinction coefficient between 0.5 and 0.7 seems to be appropriate for TMCF. In the simulation model we use the value 0.5.

Growth characteristics are assumed to be related to shade tolerance and apply for all PFTs with the same shade tolerance, except in the case of respiration, where parameters are unique to each PFT. In a Panamanian tropical lowland rain forest Ellis and others (2000) measured the maximum rate of photosynthesis (A_max) for some tree species (or species of the same genus) that also occur in the cloud forest of the study area: Trema micrantha (shade-intolerant, 18 µmol(CO₂)·m^-2s^-1), Palicourea guianensis (intermediate, 16 µmol(CO₂)·m^-2s^-1), Zanthoxylum beliziense (intermediate, 15 µmol(CO₂)·m^-2s^-1), Beilschmiedia pendula (shade-tolerant, 12 µmol(CO₂)·m^-2s^-1). Dillenburg and others (1995) report the maximum photosynthetic rates of sun leaves of Liquidambar styraciflua, which is classified as an intermediate species, to be higher than 20 µmol(CO₂)·m^-2s^-1 at an irradiance of 2000 µmol(photons)·m^-2s^-1. For the three light groups A_max is estimated to be 20 µmol(CO₂)·m^-2s^-1 for shade-intolerant PFTs, 16 µmol(CO₂)·m^-2s^-1 for intermediate and 10 µmol(CO₂)·m^-2s^-1 for shade-tolerant PFTs.

Respiration processes can be divided into growth respiration during the build-up of new biomass and maintenance respiration of living biomass. Growth respiration is estimated to amount to 25% of net production of the tree according to Ryan (1991). Parameters of maintenance respiration and the light-use efficiency (α) were fitted such that simulated maximum diameter growth for the different species (that is, under full sunlight) generated an upper envelope for measured diameter increment of TMCF species from the Botanical Garden “Francisco Javier Clavijero” in Xalapa, Mexico, and from TMCF fragments (Williams-Linera 1996; Figure S2). The transmission coefficient of leaves and the parameter for conversion of CO₂ into organic dry matter are from Larcher (2001).

For 26 species that occur in the TMCF of central Veracruz, or species of the same genus, stem wood densities are available (G. Bárcenas, unpublished data; Aguilar-Rodríguez and others 2001). Wood densities of species with the same shade tolerance were averaged. Resulting arithmetic means are 0.55 g/cm³, 0.65 g/cm³ and 0.7 g/cm³ for shade-intolerant, intermediate and shade-tolerant species, respectively. No marked differences between wood densities of the groups were found. This corresponds to the findings of Aguilar-Rodríguez and others (2001) that the majority of cloud forest species shows intermediate wood densities.

Maximum height and diameter of the different PFTs, crown length fraction and LAI of single trees were estimated (G. Williams-Linera, pers. observation). Parameters of the diameter-height relation were fitted to data from Aguilar-Rodríguez and others (2001). The form factor was adjusted such that calculated biomass corresponded to empirical biomass functions (Ter-Mikaelian and Korzukhin 1997; Acosta-Mireles and others 2002; Figure S1). The fraction of stem wood biomass was estimated to be 0.7 (compare Köhler 2000).

Basic mortality rates were defined such that pioneer species have a higher mortality compared to late successional species (for example, Poorter and Arets 2003) and that mortality decreases with increasing maximum height. In Mexican cloud forests trees die mainly through uprooting and less often through snapping or while standing (Bracho and Puig 1987; Williams-Linera 1991). In central Veracruz, the fraction of fallen dead trees was 84% (Williams-Linera 2002), in Tamaulipas 75% of the dead trees had fallen (Arriaga 1987). Most of the fallen trees belong to upper canopy species that are taller than 15 m (Arriaga 1987, 2000) or to greater diameter classes (Williams-Linera 2002). Based on these data, in the model 80% of dying trees with a dbh greater than 35 cm (PFTs 4, 5, 6) fall and damage the vegetation in the grid cell where their crown hits the ground. Remaining mortality and recruitment parameters were adjusted such that the simulated forest at steady state resembled the inventory data with respect to stem numbers, basal area and the diameter distributions of the different PFTs (Figures S3, S4).

An extensive sensitivity analysis ahs been performed with the parameterization of FORMIND for temperate rain forest in South Central Chile (Rüger and others 2007).

3. MODEL EVALUATION
Old-growth forest characteristics

To test the ability of the model to reproduce observed characteristics of old-growth forest in the study area, we used inventory data of five TMCF fragments (0.1 ha each) in central Veracruz (Williams-Linera 2002). The fragments represent the least disturbed forests in the region, but show signs of minor wood extraction and cattle grazing (Williams-Linera 2002). We grouped the individuals according to their PFT and assigned them to the corresponding diameter class (20 diameter classes of 5 cm) and used the inventory data as initial situation for the model runs. We ran 10 simulations of the dynamics of 1 ha TMCF for 400 years. In the year 400, diameter distributions for the different PFTs were compared to the inventory data. For comparison with field data, the simulation area was divided into 10 plots of 0.1 ha each and mean and standard deviation were calculated for the five forest fragments and the ten simulated plots.

Figure S3 shows the simulated development of the old-growth forest using the inventory data as initial situation. At the beginning of the simulation period, stem numbers (Figure S3A) declined due to the lack of small individuals (< 5 cm dbh) in the inventory data, before they recovered and fluctuated around the observed value. Total basal area (Figure S3B) decreased by 10 m² ha^-1 in the long term, but the relative importance of PFTs was maintained. Simulated total stem numbers and basal area were found to be in the range of observations in the forest fragments of the study area (810–1700 ind ha^-1, 35–89 m² ha^-1, Williams-Linera 2002).

Figure S4 shows stem number-diameter distributions of inventory data and simulated old-growth forest for all PFTs. For PFTs 5 and 6, which account for the largest trees and represent most of the forest’s basal area, the diameter distributions were well reproduced. Similarly, mean and variability of stem numbers of pioneer species (PFT 1) were well captured. Stem numbers for PFT 2 were underestimated, and overestimated for PFT 3 and PFT 4 in the smallest diameter class. However, these differences seem tolerable, especially when the high variability among the study sites is considered.

We validated the model by comparing simulation results regarding forest regeneration with independent field data from a chronosequence study by Muñiz-Castro and others (2006). Muñiz-Castro and others (2006) studied the arboreous vegetation in 15 abandoned pastures ranging from a few months to 80 years in age. In each site, four 10 m × 10 m (that is, 400 m²) plots at a distance of 0–10 m and 40–50 m from the border with remaining old-growth forest were sampled. Basal area, density, mean height, and mean maximum height of all individuals 5 cm or greater dbh were determined. For comparison with simulation results, we used only the data at a distance of 40–50 m from the forest border. We simulated the regeneration of TMCF starting from a treeless area of 1 ha for 100 years. We assumed that seed input is not limited, and that no further disturbances – other than gap creation by falling trees – occur during the course of succession. For comparison with field data, the simulation area was divided into 25 plots of 20 m × 20 m (that is, 400 m²). Mean and range of stem numbers, basal area, mean and maximum height for individuals 5 cm or greater dbh were calculated.

The model predicted forest regeneration reasonably well (Figure S5), taking into account that the model was calibrated such that single-tree growth and old-growth forest characteristics were reproduced. Compared to the field data, the model tends to overestimate the speed of forest recovery. During the first 10 years, differences between model simulations and field observations are due to the higher number of trees in the model. Additionally, the trees in the model seem to grow slightly faster than in reality, because maximum height is overestimated and the simulated increase of basal area is slightly faster and occurs about five years earlier than in the field data. The increase of mean height is well described by the model.

As recruitment rates were adjusted to reproduce old-growth forest conditions, they reflect regeneration after gap creation or under a closed forest canopy, and competition with non-arboreous vegetation was neglected. The data of Muñiz-Castro and others (2006) however, were sampled on abandoned pastures where native grasses dominated at the time of abandonment. This may be an explanation for the slight overestimation of stem numbers during the first years of succession. Thus, the model can be regarded to simulate an “ideal case”, that is, succession from bare ground without limitation of seed input or competition with vegetation that established during previous land use.

REFERENCES
Acosta-Mireles M, Vargas-Hernández J, Velázquez-Martínez A, Etchevers-Barra JD. 2002. Aboveground biomass estimation by means of allometric relationships in six hardwood species in Oaxaca, Mexico. Agrociencia 36: 725–736.

Aguilar-Rodríguez S, Abundiz-Bonilla L, Barajas-Morales J. 2001. Comparación de la gravedad específica y características anatómicas de la madera de dos comunidades vegetales en México. Anales del Instituto de Biología, Universidad Autónoma de México, Serie Botánica 72: 171–185.

Arriaga L. 1987. Perturbaciones naturales por la caída de árboles. Puig H, Bracho R, editors. El Bosque Mesófilo de Montaña de Tamaulipas. México, D.F.: Instituto de Ecología. p133–152.

Arriaga L. 2000. Types and causes of tree mortality in a tropical montane cloud forest of Tamaulipas, Mexico. Journal of Tropical Ecology 16: 623–636.

Bracho R, Puig H. 1987. Producción de hojarasca y fenología de ocho especies importantes del estrato arbóreo. Bracho R, Puig H, editors. El Bosque Mesófilo de Montaña de Tamaulipas. México, D.F.: Instituto de Ecología. p81–106.

Bronstein IN, Semendjajew KA. 1991. Taschenbuch der Mathematik. Stuttgart, Leipzig: Teubner.

Dillenburg LR, Teramura AH, Forseth IN, Whigham DF. 1995. Photosynthetic and biomass allocation responses of Liquidambar styraciflua (Hamamelidaceae) to vine competition. American Journal of Botany 82: 454–461.

Ditzer T. 1999. Struktur und Dynamik natürlicher und bewirtschafteter Dipterocarpaceen-wälder: Eine Fallstudie der Modellbildung und Simulation für die Bewirtschaftung natürlicher Ressourcen unter Unsicherheit. Göttinger Beiträge zur Land- und Forstwirtschaft in den Tropen und Subtropen 131. Göttingen: Erich Goltze Verlag.

Ellis AR, Hubbell SP, Potvin C. 2000. In situ field measurements of photosynthetic rates of tropical tree species: a test of the functional group hypothesis. Canadian Journal of Botany 78: 1336–1347.

Hafkenscheid R. 2000. Hydrology and biogeochemistry of tropical montane rain forests of contrasting stature in the Blue Mountains, Jamaica. PhD thesis, Vrije Universiteit Amsterdam.

Kira T. 1978. Community architecture and organic matter dynamics in tropical lowland rain forests of Southeast Asia with special reference to Pasoh Forest, West Malaysia. Tomlinson PB, Zimmermann MH, editors.Tropical Trees as Living Systems. Cambridge: Cambridge University Press. p26–30.

Köhler P. 2000. Modelling anthropogenic impacts on the growth of tropical rain forests. PhD thesis, University of Kassel, Kassel. Osnabrück: Der Andere Verlag.

Larcher W. 2001. Ökophysiologie der Pflanzen. Stuttgart: Verlag Eugen Ullmer.

Monsi M, Saeki T. 1953. Über den Lichtfaktor in den Pflanzengesellschaften und seine Bedeutung für die Stoffproduktion. Japanese Journal of Botany 14: 22–52.

Muñiz-Castro MA, Williams-Linera G, Rey-Benayas JM. 2006. Distance effect from cloud forest fragments on plant community structure in abandoned pastures in Veracruz, Mexico. Journal of Tropical Ecology 22: 431–440.

Poorter L, Arets EJMM. 2003. Light environment and tree strategies in a Bolivian tropical moist forest: an evaluation of the light partitioning hypothesis. Plant Ecology 166: 295–306.

Rüger N, Gutiérrez AG, Kissling WD, Armesto JJ, Huth A. 2007. Ecological impacts of different harvesting scenarios for temperate evergreen rain forest in southern Chile – a simulation experiment. Forest Ecology and Management 252: 52–66.

Ryan MG. 1991. Effects of climate change on plant respiration. Ecological Applications 1: 157–167.

Ter-Mikaelian MT, Korzhukin MD. 1997. Biomass equations for sixty-five North American tree species. Forest Ecology and Management 97: 1–24.

Thornley HMJ, Johnson IR. 1990. Plant and Crop Modelling – A mathematical approach to plant and crop physiology. Oxford: Clarendon Press.

Williams-Linera G. 1991. Nota sobre la estructura del estrato arbóreo del bosque mesófilo de montaña en los alrededores del campamento “El Triunfo”, Chiapas. Acta Botánica Mexicana 13: 1–7.

Williams-Linera G. 1996. Crecimiento diamétrico de árboles caducifolios y perennifolios del bosque mesófilo de montaña en los alrededores de Xalapa. Madera y Bosques 2: 53–65.

Williams-Linera G. 2002. Tree species richness complementarity, disturbance and fragmentation in a Mexican tropical montane cloud forest. Biodiversity and Conservation 11: 1825–1843.
FIGURE CAPTIONS