A spatial simulation model to explore the long-term dynamics of podocarp-tawa forest fragments, northern New Zealand

Ecological Modelling 357 (2017) 35–46

Contents lists available at ScienceDirect

Ecological Modelling journal homepage: www.elsevier.com/locate/ecolmodel

A spatial simulation model to explore the long-term dynamics of podocarp-tawa forest fragments, northern New Zealand Narkis S. Morales ∗ , George L.W. Perry School of Environment, University of Auckland, Private Bag 92019, Auckland 1001, New Zealand

a r t i c l e

i n f o

Article history: Received 21 October 2016 Received in revised form 5 April 2017 Accepted 8 April 2017 Available online 15 May 2017 Keywords: Spatially explicit individual-based model Forest fragments Long-term dynamics Long-term viability Forest management

a b s t r a c t Understanding the interactive effects of fragmentation and invasive species on forest dynamics requires a long-term perspective because they are difficult to assess in the medium- to long-term using observational or experimental data alone. In such settings ecological models have an important role to play. Here we describe the implementation of a spatially explicit individual-based model (SEIBM) representing the dynamics of small forest fragments in northern New Zealand based on empirical data collected in the region. In addition, we performed a baseline analysis to determine how well the model captured podocarp-tawa forest dynamics, and compared its performance with stand structure data obtained from an unfragmented forest in northern New Zealand. We used sensitivity analysis to determine how sensitive the model was to changes in the input parameters. In addition, we simulated different scenarios under diverse management conditions to explore the model’s potential as a management tool. The model captures the stand structural characteristics of the fragments reasonably well but under-predicts stand basal area, suggesting that it does not represent the long-term suppression of some canopy tree species adequately. Although some refinement is needed to improve its performance, we believe that the model presented here is a useful tool for management purposes and for the assessment of the long term viability of forest fragments. The model can help inform managers and decision-makers regarding the long-term persistence of podocarp-tawa forest patches. © 2017 Elsevier B.V. All rights reserved.

1. Introduction Fragmented ecosystems arise from anthropic actions such as agricultural clearing and timber extraction (Primack et al., 2001) and have resulted in the reduction of many forest ecosystems to isolated clusters of remnant habitat surrounded by a productive matrix (Santo-Silva et al., 2016; Echeverria et al., 2006). The remnant forest fragments experience changes in abiotic and biotic conditions that can have a profound influence on population and community processes such as seed dispersal, seedling establishment and ultimately species composition (Herrera and García, 2010; Broadbent et al., 2008). Understanding how human activities such as fragmentation affect forest ecosystems and how these effects might be passively or actively mitigated is challenging because the processes being considered occur over large extents in space and time and decisions are often made under high epistemic uncertainty (Messier et al., 2015). Furthermore, ecosystem

∗ Corresponding author. E-mail address: [email protected] (N.S. Morales). http://dx.doi.org/10.1016/j.ecolmodel.2017.04.007 0304-3800/© 2017 Elsevier B.V. All rights reserved.

response to the changes induced by humans (e.g. shifts in disturbance regime) may not be immediate and/or linear, but rather be lagged and/or non-linear (Johnstone et al., 2016), which makes predicting them even more challenging. In such settings it is difficult to design, implement and conduct classical experiments – rather a pluralistic approach, in which empirically grounded modelling plays important heuristic and predictive roles, is required (see Bowman et al., 2015). The podocarp-tawa forests of Northern New Zealand typify the stresses that fragmented forests suffer from and the challenges involved in understanding how to forecast and mitigate their longterm effects. This forest type is highly fragmented, under pressure from exotic herbivores and suffering invasion by exotic weeds (Ewers et al., 2006; Smale et al., 2008; Innes et al., 2010; Burns et al., 2012). Tawa (Beilschmiedia tawa; Lauraceae), a major component of the canopy of this forest type, is a long-lived (up to c. 450 y; West, 1986) endemic native tree species that is showing signs of recruitment failure (Morales, 2015; Morales et al., 2016). In response to these and other pressures, various management strategies have been applied to remnant podocarp-tawa forests. One of the most common management practices in New Zealand is the fencing of

36

N.S. Morales, G.L.W. Perry / Ecological Modelling 357 (2017) 35–46

forest fragments to exclude large herbivores, often followed by the control of small mammalian pests (Dodd et al., 2011; Burns et al., 2012). Passive activities, such as fencing, are globally the most frequently adopted restoration approach (Melo et al., 2013) and encompass a series of activities that do not require constant, active human involvement. To be successful passive approaches do, however, require the removal of activities negatively impacting the ecosystem (e.g. exclusion of large herbivores by fencing) followed by ongoing monitoring to evaluate ecosystem recovery. However, the benefits of passive restoration activities on forest fragments are difficult to assess in the medium- to long-term. In such settings ecological models can help to focus limited resources and guide management actions (Honrado et al., 2016). Ecological models are of particular value in developing and evaluating effective restoration strategies for long-lived taxa such as tree species (Anand and Desrochers, 2004). Many types of ecological models have been developed to represent vegetation dynamics across a wide range of different ecosystems (Scheller and Mladenoff, 2007; Perry and Millington, 2008). Individual-based models, or IBMs, have become one of the most widely employed types of forest model (Grimm and Railsback, 2005). An IBM represents a collection of single organisms, sometimes known as agents (e.g. a tree), their environment and the interactions between them (Reuter et al., 2011). Forest gap models, a type of IBM, were first developed in the late 1960s to study forest succession (Botkin et al., 1972; Shugart, 1984). Many existing individual-based forest gap models are derived from the JABOWA model, which was designed as a community-level forest dynamics simulator by Botkin et al. (1972). A second gap model called FORET (Shugart, 1984) was developed based on JABOWA and is the basis for 80% of gap models developed prior to 1995 (Liu and Ashton, 1995); widely used contemporary gap models derived from FORET include LINKAGES, SORTIE and FORMIND (Liu and Ashton, 1995; Deutschman et al., 1997; Dislich et al., 2009). In New Zealand, several gap model derivatives have been developed. For example, Develice (1988) developed a model called FORENZ, derived from JABOWA-FORET, but without including temperature and light limitation on growth, to predict forest dynamics on slip-faces in Fiordland. Hall and Hollinger (2000) used a model called LINKNZ, adapted from LINKAGES, to simulate conifer–hardwood and beech species forest succession. LINKNZ has also successfully been used to explore vegetation and climate shifts in the early Holocene (McGlone et al., 2011). Kunstler et al. (2009, 2013) used a model derived from SORTIE, SORTIE-NZ, to simulate the dynamics of lowland temperate rain forests on the South Island. Thrippleton et al. (2014) used a model called LANDCLIM (Schumacher et al., 2004), based on FORCLIM (Bugmann, 1994) (a descendant of the FORECE model by Kienast (1987)), to simulate long-term ecological dynamics across the large spatiotemporal scales relevant for New Zealand’s forests. However, when compared to other regions, such as North America, the use of individual-based forest models in New Zealand has been limited due to the difficulty of collecting parameterization data and the presence of unusual functional types such as long-lived pioneers (Shugart and Smith, 1996). Here we describe the implementation of a spatially explicit individually-based model (SEIBM) designed to represent the longterm dynamics and viability of podocarp-tawa forest fragments in New Zealand. Our main purpose was to design and implement a model that could simulate the fate of forest fragments across a range of management contexts, representing the effects of both fragmentation and fencing. Also, as a proof-of-concept we simulated different scenarios under diverse management conditions to explore the model’s potential usefulness as a management tool.

2. Methods 2.1. Model description The model described in this section is a spatially explicit individually-based model (SEIBM) designed to represent the longterm dynamics and viability of podocarp-tawa forest fragments under different types of management (unfragmented forest, fenced and unfenced fragments) of a forest of 16 ha (i.e. a 100 × 100 grid of 4 × 4 m cells) (Fig. 1). The model design is, by intention, reasonably generic and could be adapted to model other types of forests if needed. The model was essentially manually calibrated but wherever possible we drew on empirical estimates in its parameterization and we evaluated the model’s outcomes against forest stand structures described in the NZ literature (although such information was not available for all parameters or outcomes). The information provided in this section is a summary of the Overview, Design concepts, and Details (ODD) protocol described by Grimm et al. (2006, 2010) as presented in Appendix A in Supplementary material 1. The model was implemented using NetLogo 5.3.1 (Wilensky, 1999); a free and open source software tool, available under a GPL license. The model considers three forest canopy species – Beilschmiedia tawa, Laurelia novae-zelandiae and Dacrydium cupressinum – and three sub-canopy species that typify these forest ecosystems: Hedycarya arborea, Piper excelsum and Melicytus ramiflorus. In each time-step or tick (one year) the model sequentially processes a series of ecological routines as follows: restoration planting, seed dispersal (both from within and beyond the fragment), herbivory, mortality, gap formation, regeneration (recruitment), growth and seedling/sapling demography (Fig. 1) The model was initialized with a spatially random pattern (there are not data describing individual-level spatial pattern in these forests) and parameterized using data collected from unfragmented podocarp-tawa forest (continuous fenced forest) and from previously published literature. The study area for data collection was located near Cambridge in the Waikato region, North Island, New Zealand. The unfragmented forests were located at Te Miro Scenic Reserve, Maungatautari Ecological Island North and Maungatautari Ecological Island South (Morales et al., 2016). The Maungatautari Ecological Island is surrounded by a predator-proof fence (Burns et al., 2012) (Fig. 2). At initialization the species abundances, and the size (diameter at breast height; dbh) and the age of each individual adult tree will depend on the scenario being considered (unfragmented forest, fenced and unfenced fragment). Regardless of the scenario, the model starts with the same number of seedlings of each species (10 seedlings [individuals < approximately 50 cm in height]) regardless of the scenario and saplings (one sapling [individuals > approximately 50 cm in height) in each grid cell. When simulated, restoration activities (i.e. planting) are represented by adding a species-specific amount of saplings to the sapling bank at a specific frequency. Each time-step, adult trees produce a number of seedlings following a Poisson deviate specified by a seed production parameter (seed-prod). We distinguish between three types of dispersal (Fig. 3): under individual trees (‘neighbourhood dispersal’), within the fragment (‘local dispersal’) and from beyond the fragment (‘long-distance dispersal’; by, for example, frugivorous birds). Neighbourhood dispersal represents the dispersal that occurs under parent trees. The area within which this occurs depends on the size of the individual tree’s crown. Local dispersal represents those seedlings that establish within the fragment but at some distance from their parent. The process of long-distance seed dispersal is intended to give a chance for all the species to invade and establish a patch during the succession via movement of seeds

N.S. Morales, G.L.W. Perry / Ecological Modelling 357 (2017) 35–46

37

Fig. 1. Flow-chart representing the sequence of the ecological processes represented in the SEIBM; dashed lines represent processes that can be turned on and off.

from beyond the fragment. The proportion of seeds actively dispersed by each species depends on how likely their seeds are to be dispersed by frugivorous birds (a number of NZ’s late successional forest species, including B. tawa, are bird-dispersed). New seedlings are added to each grid cell’s seedling bank each year. Saplings can

appear in a grid cell either as seedlings transitioning to the sapling stage or via the restoration planting process. Annual tree (adult) mortality follows the formula of Shugart (1984), which yields a background mortality rate such that 98% of individuals die before reaching the species maximum age. In

38

N.S. Morales, G.L.W. Perry / Ecological Modelling 357 (2017) 35–46

Fig. 2. Map of the study site, near Cambridge in the Waikato region, North Island, New Zealand. Thick black lines represent roads. Polygons in light green are forest fragments in the study area. Te Miro, Maungatautari north and Maungatautari south, in dark green, are areas where the data were collected to parameterize the model. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article).

Fig. 3. Schematic description of the dispersal processes: (a) neighbourhood dispersal, (b) local (within fragment) dispersal and (c) long-distance dispersal (from beyond the fragment).

N.S. Morales, G.L.W. Perry / Ecological Modelling 357 (2017) 35–46

39

Fig. 4. Examples of different forest fragments, (a) example of a fenced fragment (b) example of an unfenced fragment and (c) photo of a forest fragment boundary. All fragments are located in Te Miro, Waikato region, North Island, New Zealand.

the absence of data describing age-dependent mortality rates for all species and size-classes that we simulate (Richardson et al., 2009 present rates for B. tawa and D. cupressinum, but only for large individuals) we adopted the most parsimonious model of ageindependent background mortality. If the dying individual belongs to a gap-maker species then on its death it creates a canopy gap, potentially extending beyond its own grid cell and affecting its neighbours. Each species growth is represented by an annual diameter increment (D; m/yr), calculated using the formula provided by Botkin et al. (1972) for the JABOWA gap model. The height (in meters) of a tree of known dbh is calculated using the formula provided by Ker and Smith (1955), which is the standard approach for gap models and their derivatives. Diameter increment is corrected by a neighbourhood effect that represents competition and an edge effect. The neighbour effect (competitive penalty), phenomenologically represents the intensity of competition that each individual is experiencing and is based on its height relative to those of its neighbours, which is a proxy for local light availability (our approach is based on that of Dislich et al., 2009).

2.2. Baseline analysis We explored the dynamics of an unfragmented forest of 16 ha (i.e. a 100 × 100 cell grid of 4 × 4 m cells) by running the model 30 times for 1000 years under baseline parameter conditions. We graphically analysed three of the model’s key state variables: the density of the dominant late-successional species (B. tawa, L. novaezelandiae and D. cupressinum) included in the model, their mean dbh and mean basal area. We did not use inferential statistics to evaluate the model’s dynamics and behaviours. Simulation-based analyses can be highly replicated, which leads to very small p values regardless of the effect size and, ultimately, ‘guaranteed’ rejection of null hypotheses of no difference between treatments (White et al., 2014).

2.3. Sensitivity analysis Given the uncertainties in the parameter estimations we conducted a sensitivity analysis to assess how sensitive the model is to its parameterization and hence assess the plausibility of the model’s dynamics and the confidence we can have in its outcomes. For the sensitivity analysis we used 12 input parameters and three state variables. The parameters assessed for their sensitivity were external-species, repro-age, sapling-survival, seedling-survival, seedling-survival, ldd-dispersal-dist, ldd-dispersalfrac, regen-height, seed-prod, shade-tolerance, supp-mortality and sup-tolerance (Appendix B in Supplementary material 2; Table A2.1). We chose these parameters because we deemed them a priori to be the most important ecological traits and many of them are not well quantified in the literature. The outputs (state variables) we considered in the sensitivity analysis were proportional abundance, mean dbh, and mean basal area (m2 /ha) of all species (Appendix B in Supplementary material 2; Table A2.1). Based on the point at which the composition of the simulated forest stabilises, 1000 years was deemed an appropriate simulation length for the sensitivity analyses (Fig. 5). However, D. cupressinum did not stabilise after 1000 years (especially in terms of density and mean basal area); rather the simulated forest was characterised by a very few, very large individuals of this species and so its abundance and mean size vary continually. This pattern is consistent with the known dynamics of this species, which can be seen as an extremely long-lived pioneer (Ogden and Stewart 1995). We changed the baseline values by ±20% for all species, one parameter at a time, and ran the model 30 times for 1000 years for all 12 parameters. Thus, the sensitivity analysis constitutes a local ‘oneat-a-time’ approach and considers neither parameter interactions nor spatial pattern (Hamby, 1994; Drechsler, 1998). To estimate each parameter’s sensitivity, the mean for each state variable across the 30 replicate runs of the model for each param-

40

N.S. Morales, G.L.W. Perry / Ecological Modelling 357 (2017) 35–46

fragment and an unfenced fragment both with and without restoration.). We used a timeframe of 500 years in a forest fragment of 16 ha (i.e. a 100 × 100 grid of 4 × 4 m cells). We used the same area as in our baseline analysis for ease of comparison of the changes in the simulated fragment, but in reality most podocarp-tawa fragments are considerably smaller than 16 ha. We used field data described in Morales et al. (2016) to initialize the different scenarios. To simulate the dynamics of a fenced fragment only the edge effect was ‘on’ (b2 parameter set to 0.15 representing an approximately 30 m edge effect, see ODD description in SM). The unfenced fragment was a combination of edge effect and herbivory routines. Restoration planting was simulated by adding one sapling per grid cell of B. tawa, L. novae-zelandiae and D. cupressinum to the sapling bank every five years; this equates to a planting density of 100 individuals per species per hectare. Finally, dispersal from beyond the fragment was deactivated for all the scenarios to simulate a completely isolated forest fragment. Fig. 5. Species density after 1000 years across n = 30 model realizations. Coloured lines represent densities of species using field data obtained from an unfragmented forest (red: B. tawa, green: L. novae-zelandiae and orange: D. cupressinum). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article).

eter of interest was calculated, and sensitivity quantified following Hamby (1994): S y,x = (y/yb ) (x/xb )

(9)

where: Sy,x is the sensitivity of the state variable y with respect to change in parameter x, x and y are the absolute change in x and y between the baseline and the change in baseline values (±20%) and xb and yb are the values under baseline conditions. We deemed state variables (y) with Sy,x > 1 as sensitive to changes in parameter (x) (Hamby, 1994). For each species we ranked the parameters’ sensitivity values (Sy,x ) according to the response of the state variables (Appendix B in Supplementary material 2) and we then averaged the rankings of each species to generate an overall rank sensitivity for each parameter (mean of the mean rank). 2.4. Scenario-based evaluation of the long-term dynamics We evaluated two management activities: fencing and restoration activities (e.g. seedling planting) in a factorial way (i.e. a fenced

3. Results 3.1. Baseline analysis Across the six species, adult tree density is similar across replicate simulations; the density of B. tawa and L. novae-zelandiae differs between runs, but only slightly (around 5%). This dynamic is unsurprising because B. tawa and L. novae-zelandiae have similar habitat requirements so they are constantly competing for resources and tend to dominate late-successional (low light) environments (Figs. 5 and 6). The spatial dynamics of these two species suggest a localised lock-in; when a species starts to do well it will capture more grid cells, increase in density, minimise interspecific neighbourhood interactions, so capture more cells, and so forth. We compared the densities of B. tawa, L. novae-zelandiae and D. cupressinum after 1000 simulation years with the field estimates for unfragmented forest described in Morales et al., 2016. For B. tawa the model yielded a mean density of 291 individuals per ha (95% CI [286,295]) while the field data yielded 154 individuals per ha (95% CI [135,176]). L. novae-zelandiae mean density were 293 individuals per ha (95% CI [288,298]) and 134 individuals per ha (95% CI [118,154]) for the model and field data, respectively. For D. cupressinum the model gave a mean density of 22 individuals per ha (95% CI [21,23]) while the field data estimates two individuals per ha (95% CI [2,3]).

Fig. 6. Species densities (stems/ha) in (a) from n = 10 randomly selected model realizations under baseline conditions, each of 1000 years. The black box represents the point where the densities of the dominant late-successional species stabilized (based on visual assessment). The map in (b) shows the composition and spatial configuration of an unfragmented area of forest (400 × 400 m; 16 ha) after 1000 years of simulation. Black represents gaps in the forest (i.e. cells without an adult or juvenile individual). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article).

N.S. Morales, G.L.W. Perry / Ecological Modelling 357 (2017) 35–46

Fig. 7. Species mean dbh (m) in (a) for each of the six species in n = 10 randomly selected model realizations under baseline conditions (same realizations from Fig. 6), across 1000 years; colouring as per Fig. 6 . The erratic trajectories for D. cupressinum occur because at any time there are just a few large individuals that produce an abrupt drop in the species mean dbh when they die. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article).

41

Fig. 8. Species mean basal area (m2 /ha) after 1000 years (n = 30 model realizations). Coloured horizontal lines represent mean basal area of species for B. tawa, L. novae-zelandiae and D. cupressinum estimated from field data obtained from the unfragmented forest. Colour as per Fig. 5. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article).

Table 1 Comparison of the percentage of individuals per size-class given by the model (mean ± SE; n = 30) and the observed data for B. tawa and L. novae-zelandiae. Model data are based on n = 30 model realizations and observed data are based on n = 62 for B. tawa and n = 54 for L. novae-zelandiae. B. tawa Model

L. novae-zelandiae Observed

Model

Observed

Size-class (dbh) Individuals (%) Individuals (%) Individuals (%) Individuals (%) 0–20 cm 20–40 cm 40–60 cm 60–80 cm 80–100 cm 100–120 cm >120 cm

46.8 ± 4.8 25.9 ± 3.1 13.4 ± 1.9 7.9 ± 1.4 4.2 ± 1.4 1.9 ± 0.6 0

35 29 19 5 10 0 2

45.8 ± 4.5 26.2 ± 3.2 13.9 ± 1.9 8.0 ± 1.2 4.3 ± 0.8 1.8 ± 0.6 0

43 30 19 4 4 0 2

Mean dbh for both B. tawa and L. novae-zelandiae are stable during each model run and show little variation between simulations, a dynamic consistent with the fact that density trajectories vary little between model runs (Fig. 7). For B. tawa the model estimated a mean dhb of 23 cm (95% CI [22.8, 23.2]) while we obtained a dbh of 36.1 cm (95% CI [29.4, 42.6]) from the field data (Table 1). For L. novae-zelandiae the model estimated a mean dhb of 23 cm (95% CI [22.9, 23.2]) while a dbh of 31.1 cm (95% CI [23.8, 38.3]) was obtained from the field data (Table 1). For D. cupressinum the model gave an estimated mean dbh of 38 cm (95% CI [37.4, 38.9]). We could not compare the model results for D. cupressinum with the field data obtained from the unfragmented forests as we recorded only one individual (with a diameter of 1.7 m). However, previous studies have reported that dbh for D. cupressinum varies considerably from site-to-site with most of the individuals in the range of 30–45 cm (Stewart et al., 1998; Rogers, 1999 – both South Island studies). We compared the mean basal area results from the model after 1000 years with those estimated in the field (Figs. 8 and 9). For B. tawa the model estimated a stand-level mean basal area of 26.5 m2 /ha (95% CI [25.9, 27]) while the field data results yielded 24 m2 /ha (95% CI [21.1, 27.5]). For L. novae-zelandiae the estimated mean basal area was 25.8 m2 /ha (95% CI [25.2, 26.5]) while the field data yielded 17.9 m2 /ha (95% CI [15.6, 20.4]). D. cupressinum estimated mean basal area was 5.9 m2 /ha (95% CI [5.6, 6.2]) as compared to 6.1 m2 /ha (95% CI [5.3, 6.9]) for the field data.

Fig. 9. Species mean basal area in for each of the species across n = 10 randomly selected realizations (same realizations as in Fig. 6), over 2000 years; colouring as per Fig. 6. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article).

3.2. Sensitivity analysis The angiosperm canopy species B. tawa and L. novae-zelandiae were insensitive to changes in the 12 parameters that we analyzed. However, the gymnosperm canopy tree species D. cupressinum was sensitive to sapling survival, seedling survival and regeneration height parameters. For the sub-canopy species, the most sensitive parameters for H. arborea were seedling survival, fraction of dispersed seedlings and suppression mortality (see Appendix A in Supplementary material 1; section 3.3.5). M. ramiflorus was sensitive to sapling survival, seedling survival and suppression mortality (Table 2). P. excelsum was not included in Table 2 as one of the state variables (mean basal area) showed no changes in the sensitivity index, and so was not comparable with the other species in that table. Proportional density and mean dbh were most sensitive to changes in the fraction of seedlings dispersed within the fragment, regeneration height and seed production. Overall, seedling survival was the parameter with the highest sensitivity rank across all the analysed species (excluding P. excelsum; Table 2). The complete results of the sensitivity ranking exercise are summarized in Appendix B in Supplementary material 2.

42

N.S. Morales, G.L.W. Perry / Ecological Modelling 357 (2017) 35–46

Table 2 Overall ranks of the parameters’ sensitivity to a change of ± 20% in all species parameters simultaneously. The parameters are described in Appendix B in Supplementary material 2. Species are (a) H. arborea, (b) M. ramiflorus and (c) D. cupressinum. Note that B. tawa and L. novae-zelandiae were insensitive to all 12 input parameter so are not included here. Species

Rank of the mean rank

Parameter

(a)

(b)

(c)

seedling-survival sapling-survival ldd-dispersal-frac seed-prod supp-mortality regen-height repro-age seedling-survival ldd-dispersal-dist shade-tolerance supp-tolerance external-species

1 9 3 4 2 11 11 4 7 7 4 9

3 1 6 4 2 6 4 6 6 10 12 10

1 2 4 8 11 4 3 9 4 7 9 11

1 2= 2= 4= 4= 6 7= 7= 7= 7= 11 11

B. tawa mean dhh (m) was 21.5% lower than the L. novaezelandiae. In terms of basal area (m2 /ha), B. tawa had a 77% higher basal area than L. novae-zelandiae. The high Mean dhh (m), basal area (m2 /ha) and low density of D. cupressinum suggest the persistence of very few large individuals. Sub-canopy species, H. arborea and M. ramiflorus small dbh and lower basal areas for both species (Figs. 10c and 11b). In the unfenced scenario with restoration activities, B. tawa and L. novae-zelandiae dhh were almost identical. However, B. tawa had a 28.5% higher basal area than L. novaezelandia. In the case of D. cupressinum the mean dhh (m) did not vary although the basal area (m2 /ha) was 32% lower than in the previous scenario (Figs. 10d and 11c). The spatial pattern of the unfenced fragment showed that the fragment was dominated mostly by the canopy dominant species B. tawa. There are few individuals of D. cupressinum and H. arborea distributed throughout the forest fragment. M. ramiflorus tends to dominate the edges of the forest fragments (Fig. 12c). In the scenario of a fenced fragment with restoration activities seedling planting changed the composition and spatial distribution dominated only by the canopy species B. tawa and L. novae-zelandiae (Fig. 12d).

3.3. Scenario-based evaluation of the long-term dynamics of the fragments At the end of the 500-year runs B. tawa and L. novae-zelandiae were the dominant species in the fenced forest fragment with high densities. D. cupressinum and M. ramiflorus were still present at the end of the modelled period although in low densities. H. arborea and P. excelsum were short-lived and began to decline with a few decades (Figs. 10a and 11a ). When restoration activities are introduced (seedling planting), B. tawa and L. novae-zelandiae become the dominant species as in the fenced forest fragment. However, in this scenario D. cupressinum density had a substantial increase (4.5 times) compared to the scenario without restoration. The remainder of the species had an abundance of less than one individual per ha (Figs. 10b and 11a). Mean dhh (m) and basal areas (m2 /ha) of B. tawa and L. novaezelandiae are similar, suggesting that individuals of these species are of similar size. The high mean dhh (m), basal area (m2 /ha) and low density of D. cupressinum suggest the persistence of scattered, large, individuals. Despite the higher density in comparison with D. cupressinum, M. ramiflorus have a smaller mean dbh (m) and basal area (m2 /ha), meaning that individuals of that species are small (Figs. 10a and 11b). In the fenced scenario with restoration activities basal areas (m2 /ha) of B. tawa and L. novae-zelandiae decreased 8% and 22% relative to the fenced fragments. D. cupressinum basal area increased by 51% when restoration activities occurred; this effect is due to an increase in light due to the edge effect that improves D. cupressinum growth (Figs. 10c and 11c). The spatial pattern of the fenced fragment after 500 years showed a strong edge effect with the establishment of canopy dominant species inhibited (B. tawa and L. novae-zelandiae). M. ramiflorus tends to dominate the edges of the forest fragments, while H. arborea is present in the forest fragment border and scattered around the simulated fragment (Fig. 12a). In the fenced fragment with restoration activities scenario, seedling planting changed the composition such that it is dominated by the large canopy species (Fig. 12b). At the end of the 500-year simulation of the unfenced fragment scenario B. tawa was the main dominant species with a density 7 times higher than L. novae-zelandiae. D. cupressinum, H. arborea and M. ramiflorus were highly underrepresented (Fig. 10, unfenced fragment; Fig. 11c). When restoration activities were ‘on’, B. tawa and L. novae-zelandiae became the dominant species. D. cupressinum density did not vary in comparison with the previous scenario. None of the sub-canopy species were present at the end of the 500 years simulation period (Figs. 10d and 11a).

4. Discussion 4.1. Plausibility of the model’s baseline dynamics We have described the implementation of a spatially explicit individually-based model (SEIBM) designed to represent the longterm dynamics and viability of podocarp-tawa forest in New Zealand. Our main purpose was to outline how the model represents fundamental ecological processes and how decisions about that representation were made. We performed a baseline analysis to assess how well the model captured the dynamics of podocarp-tawa forest and compared its performance with stand structure data obtained from an unfragmented forest in northern New Zealand and derived from the literature, alongside a sensitivity analysis intended to determine how sensitive the model was to changes in the input parameters. Finally, we simulated forest change under a range of management conditions to explore the model’s potential as a management tool There have been few attempts to develop individual-based forest models to represent the dynamics of New Zealand’s indigenous forests; FORENZ, SORTIE/NZ and LINKNZ (Develice, 1988; Kunstler et al., 2009; Hall and Hollinger, 2000) represent the only three such efforts to date. These models were mostly developed for forests where Nothofagaceae (Fuscaspora spp. and Lophozonia sp.) are the dominant taxa (FORENZ and SORTIE/NZ). LINKNZ provides the most comprehensive effort to test the applicability of a single gap-model across a wide range of species, and has provided insight into New Zealand conifer–hardwood and beech (Nothofagaceae) forest succession. Our model is one of the first attempts to simulate the dynamics of podocarp-tawa forests and we explicitly designed it as a tool to inform long-term conservation strategies and restoration management activities. As such, we aimed to make the model simple enough to be modified to simulate the dynamics of other forest types and to be able to be parameterized from widely available or easily collectable data. In the following discussion, we focus on the performance of the model in the context of the available field-data for the forest type we consider. In general, our model’s predictions were plausible when assessed against field data. To parameterize our model we used a broad range of the available information regarding New Zealand forest dynamics, and podocarp-tawa forests in particular; most of the parameters could be estimated reasonably accurately. We endeavoured to incorporate all available information and used

N.S. Morales, G.L.W. Perry / Ecological Modelling 357 (2017) 35–46

43

Fig. 10. Long-term simulation (500 years) of a (a) fenced without and (b) with restoration, and a (c) unfenced forest fragment without and (d) with restoration. The graphs represent the number of adult trees per ha, mean dbh (m) and mean basal area (m2 /ha) in the fenced and unfenced forest fragment. Colouring as per Fig. 6. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article).

Fig. 11. Long-term simulation (500 years) of a fenced without and with restoration, and an unfenced forest fragment without and with restoration. The bar plots represent the (a) total number of adult trees per ha ± SD, (b) total mean dbh (cm ± SD) and (c) total mean basal area (m2 /ha ± SD) in the fenced and unfenced forest fragment.

44

N.S. Morales, G.L.W. Perry / Ecological Modelling 357 (2017) 35–46

Fig. 12. The different maps shows the composition and spatial configuration of a forest block after 500 years of simulation (400 × 400 m; 16 ha). (a) represents a fenced forest fragment, (b) represents a fenced forest fragment with restoration, (c) represents an unfenced forest fragment and (d) represents an unfenced forest fragment with restoration Black cells represent gaps in the forest (i.e. cells without an adult or juveniles individual). Colouring as per Fig. 6. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article).

sensitivity analyses to confirm that the parameterisation was reasonable as were the model’s dynamics. West (1986) showed that B. tawa has higher canopy cover and tends to be the dominant species in unfragmented podocarp-tawa forest, and, conversely, that D. cupressinum has a lower canopy cover and is less abundant in the unfragmented forest. Initially, B. tawa and L. novae-zelandiae are randomly distributed across the modelled forest but they slowly shift to a more aggregated pattern. D. cupressinum develops a scattered pattern comprising a few, large emergent individuals, a dynamic characterized in the ‘lozenge model’ that describes successional change in these forest ecosystems (Ogden and Stewart 1995) (Fig. 6b). When compared to data collected from unfragmented forests the model predicts the proportional density across the different tree species reasonably well. Although the model slightly under-predicted the density of medium size-class individuals (e.g. 20–60 cm) and over-predicted the small size class, especially for B. tawa (e.g. 0–20 cm), other size classes were adequately predicted, set against another study with a higher number of samples (Table 1 cf West 1986).

The discrepancy in densities between what was obtained in the model and the field data gathered from the field data probably arises from three main factors. First, not all the canopy species present in the field were represented in the model, although those that were not (e.g. Knightia excelsa, Litsea calicaris, Rhopalostylis sapida, Alectryon excelsus) were typically present only in low densities (representing 11.5% of the total individuals per ha). Second, the overall stem density in the forest is difficult to represent on a fixed grid, as in the field the distance between individuals is not constant. In other words, the distance between individuals could be smaller or larger than the imposed grid cell size in the model (grid cells are 4 × 4 m and only one adult tree can occupy any given grid cell). For example, the model can simulate a maximum of 625 individuals per ha compared to the 418 individuals per ha from the field data (only considering the six modelled species). If we consider all of the species in the field data the density increases to 540 individuals per ha, which is still one-third lower than the model’s maximum density. Our field data show that the mean distance between individuals was 4.1 ± 2.4 m. Thus, although the mean distance between individuals is similar to the grid resolution we use, it is highly vari-

N.S. Morales, G.L.W. Perry / Ecological Modelling 357 (2017) 35–46

able. Finally, the forests we have collected empirical data from are not 1000 years old nor are they entirely undisturbed; rather they have experienced forms of anthropic and natural disturbances not included in the model. On the whole, nevertheless, the proportional densities of the different species in the model are consistent with those observed in the field (Fig. 5; Table 1). The mean tree diameters produced by the model are consistently lower than the reference unfragmented forest. However, the model’s results do not seem implausible as mean dbh from the field data is variable and context-dependent; for example, lower mean dbh values for B. tawa of 26.3 ± 2.3 cm have been reported in other B.tawa forests in northern New Zealand (Smale et al., 2008). The model dynamics do suggest, however, that simulated mortality rates are too high for smaller individuals and this might be improved by adding an age-dependent mortality sub-model (if sufficient data were available). The mean basal area of L. novaezelandiae was over-estimated in comparison with the observed data and the basal area values presented in other studies (West 1986). Overall, our model, does however, manage to represent the key structural and successional patterns observed in the reference unfragmented forest and has the potential for future improvements in the areas highlighted above (especially in the representation of suppression and mortality dynamics). The discrepancy between model predictions and field observations is, at least in part, due to the lack of regeneration in B. tawa at some sites as reported by Burns et al. (2011), Morales (2015), Morales et al. (2016). In addition, the discrepancies between model predictions and the available empirical data could arise from differences in the site characteristics (e.g. soil type, past history), which were not represented in the model. In addition, because New Zealand’s native trees, including tawa, can tolerate extremely long (potentially centennial) suppression periods (Smale et al., 1986 and Hall and Hollinger, 2000) they are challenging to represent in forest simulators. Therefore, we likely underestimated suppression tolerance, especially for smaller individuals. 4.2. Sensitivity analysis Our sensitivity analysis demonstrated that the most sensitive parameters are, in general, those for which there was the least empirical information. While there has been increased use of sophisticated reverse inference parameterization for forest simulation models, (Hartig et al., 2012, 2014) such approaches are computer intensive and require more data than we have available. Parameter estimation for seedling survival, seedling transition and sapling survival was extremely difficult as there is little, if any, reliable species-specific information available for these. Nevertheless, on the basis of the model evaluations we conducted we are confident that if there was some mis-estimation of those parameters it does not unduly affect the robustness of the model’s outcomes.

45

became the only dominant canopy species. Although restoration activities seem to help in terms of composition, the fragments are totally dependent on the influx of new seedlings from planting. The most positive results with regards to the long-term survivorship and persistence of the canopy tree species were obtained when fencing was coupled with an absence of herbivory and restoration activities such as planting were in place, especially for D. cupressinum

4.4. Model refinements As more empirical information becomes available the model could be structurally refined and its parameterization improved. For example, age-dependant mortality could improve the model substantially, especially for those species capable of persisting over long periods in low light conditions; unfortunately there is a lack of quantitative information regarding this process for the studied species other than, to some extent, B. tawa (Richardson et al., 2009; Smale et al., 1986). Furthermore, the lottery mechanism needs to be improved because stand dynamics under complete regeneration failure are inadequately represented. In those scenarios where there is high herbivory pressure and a severe edge effect some simulations can drift towards a monodominant forest community because by chance alone a sapling will appear in a cell and make the transition to adulthood; this will eventually drive a lock-in effect between a species local density and its regeneration success. Finally, the addition of a module that accounts for individual tree’s response to soil or climate conditions would allow future researchers to determine the potential effect of fertiliser spill-over and climate change on the fragments we consider. Currently there is a lack of data to achieve this. Although Smale et al. (2014) provide growth data for B. tawa, it is highly-biased in terms of specific soils. Furthermore, including ways to export the data to help better visualization of the results will be helpful. For example, combining the results with a 3D visualization tool as demonstrated by Dislich et al. (2009). The use of NetLogo as an implementation platform means that a range of user-friendly graphical tools are available to conduct model-based experiments. However, despite the plasticity of the model in some circumstances changes to the model will be needed. For example, in cases where there is more information for some of the modelled processes, such as mortality or growth, a change in the code will be needed. Although this processes could be fairly trivial for a modeller, it could quite daunting for the end user which can impact on the adaptability of the model to others forest types.

5. Conclusions 4.3. Scenario assessment via in silico experiment The outcomes of the in silico experiments suggest that fencing is necessary to ensure the long-term persistence of the dominant tree species (B. tawa, L. novae-zelandiae and D. cupressinum) in forest fragments because it greatly reduces herbivore pressure from large mammals (e.g. cattle and sheep). However, the edge effect has an impact on the composition and structure of the forest, especially in the canopy species that cannot colonize the fragment edges. The outcome of this effect is the smaller dbh and lower basal areas in smaller fragments compared with the baseline simulations (unfragmented forest). In the case of unfenced fragments, the fragments shifted towards an almost mono-dominant forest, which we deemed to represent fragment collapse. This collapse started after around 130 years at which point there was a shift and B. tawa

Our spatially explicit individually-based model (SEIBM) is capable of representing the long-term dynamics and viability of podocarp-tawa forest in New Zealand. Although the model can be improved (as any model can), we believe it is a useful tool and one capable of accomplishing the objectives for which it was designed. Areas where improvement is most needed are the representation of mortality and suppression and forest dynamics under complete regeneration failure. Nevertheless, preliminary scenario simulations show that long-term survivorship and persistence of the canopy tree species can be achieved when fencing, total eradication of herbivores and restoration activities (e.g. planting) are in place. Our results show that models can be a powerful tool than can inform management activities and the assessment of the long-term viability of forest fragments.

46

N.S. Morales, G.L.W. Perry / Ecological Modelling 357 (2017) 35–46

Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.ecolmodel.2017. 04.007. References Anand, M., Desrochers, R.E., 2004. Quantification of restoration success using complex systems concepts and models. Restor. Ecol. 12, 117–123. Botkin, D.B., Janak, J.F., Wallis, J.R., 1972. Some ecological consequences of a computer model of forest growth. J. Ecol. 60, 849–872. Bowman, D.M.J.S., Perry, G.L.W., Marston, J.B., 2015. Feedbacks and landscape-level vegetation dynamics. Trends Ecol. Evol. 30, 255–260. Broadbent, E.N., Asner, G.P., Keller, M., Knapp, D.E., Oliveira, P.J.C., Silva, J.N., 2008. Forest fragmentation and edge effects from deforestation and selective logging in the Brazilian Amazon. Biol. Conserv. 141, 1745–1757. Bugmann, H., 1994. On the Ecology of Mountainous Forests in a Changing Climate. Swiss Federal Institute of Technology, Zurich, Switzerland, Ph.D. thesis. Burns, B.R., Floyd, C.G., Smale, M.C., Arnold, G.C., 2011. Effects of forest fragment management on vegetation condition and maintenance of canopy composition in a New Zealand pastoral landscape. Austral Ecol. 36, 153–166. Burns, B., Innes, J., Day, T., 2012. The use and potential of pest-proof fencing for ecosystem restoration and fauna conservation in New Zealand. In: Somers, M.J., Hayward, M. (Eds.), Fencing for Conservation. Springer, New York, pp. 65–90. Deutschman, D.H., Levin, S.A., Devine, C., Buttel, L.A., 1997. Scaling from trees to forests: analysis of a complex simulation Model. Science 277, 1684. Develice, R.L., 1988. Test of a forest dynamics simulator in New Zealand. N. Z. J. Bot. 26, 387–392. Dislich, C., Günter, S., Homeier, J., Schröder, B., Huth, A., 2009. Simulating forest dynamics of a tropical montane forest in South Ecuador. Erdkunde 63, 347–364. Dodd, M., Barker, G., Burns, B., Didham, R., Innes, J., King, C., Smale, M., Watts, C., 2011. Resilience of New Zealand indigenous forest fragments to impacts of livestock and pest mammals. N. Z. J. Ecol. 35, 83–95. Drechsler, M., 1998. Sensitivity analysis of complex models. Biol. Conserv. 86, 401–412. Echeverria, C., Coomes, D., Salas, J., Rey-Benayas, J.M., Lara, A., Newton, A., 2006. Rapid deforestation and fragmentation of Chilean temperate forests. Biol. Conserv. 130, 481–494. Ewers, R.M., Kliskey, A.D., Walker, S., Rutledge, D., Harding, J.S., Didham, R.K., 2006. Past and future trajectories of forest loss in New Zealand. Biol. Conserv. 133, 312–325. Grimm, V., Railsback, S.F., 2005. Individual-Based Modeling and Ecology. Princeton University Press, Princeton. Grimm, V., Berger, U., Bastiansen, F., Eliassen, S., Ginot, V., Giske, J., Goss-Custard, J., Grand, T., Heinz, S.K., Huse, G., Huth, A., Jepsen, J.U., Jørgensen, C., Mooij, W.M., Müller, B., Pe’er, G., Piou, C., Railsback, S.F., Robbins, A.M., Robbins, M.M., Rossmanith, E., Rüger, N., Strand, E., Souissi, S., Stillman, R.A., Vabø, R., Visser, U., DeAngelis, D.L., 2006. A standard protocol for describing individual-based and agent-based models. Ecol. Model. 198, 115–126. Grimm, V., Berger, U., DeAngelis, D.L., Polhill, J.G., Giske, J., Railsback, S.F., 2010. The ODD protocol: a review and first update. Ecol. Model. 221, 2760–2768. Hall, G.M.J., Hollinger, D.Y., 2000. Simulating New Zealand forest dynamics with a generalized temperate forest gap model. Ecol. Appl. 10, 115–130. Hamby, D.M., 1994. A review of techniques for parameter sensitivity analysis of environmental models. Environ. Monit. Assess. 32, 135–154. Hartig, F., Dyke, J., Hickler, T., Higgins, S.I., O’Hara, R.B., Scheiter, S., Huth, A., 2012. Connecting dynamic vegetation models to data – an inverse perspective. J. Biogeogr. 39, 2240–2252. Hartig, F., Dislich, C., Wiegand, T., Huth, A., 2014. Technical Note: approximate Bayesian parameterization of a process-based tropical forest model. Biogeosciences 11, 1261–1272. Herrera, J.M., García, D., 2010. Effects of forest fragmentation on seed dispersal and seedling establishment in ornithochorous trees. Conserv. Biol. 24, 1089–1098. Honrado, J.P., Pereira, H.M., Guisan, A., 2016. Fostering integration between biodiversity monitoring and modelling. J. Appl. Ecol. 53, 1299–1304. Innes, J., King, C.M., Bridgman, L., Fitzgerald, N., Arnold, G., Cox, N., 2010. Effect of grazing on ship rat density in forest fragments of lowland Waikato, New Zealand.N. Z. J. Ecol. 34, 227–232. Johnstone, J.F., Allen, C.D., Franklin, J.F., Frelich, L.E., Harvey, B.J., Higuera, P.E., Mack, M.C., Meentemeyer, R.K., Metz, M.R., Perry, G.L., Schoennagel, T., Turner, M.G., 2016. Changing disturbance regimes, ecological memory, and forest resilience. Front. Ecol. Environ. 14, 369–378. Ker, J.W., Smith, J.H.G., 1955. Advantages of the parabolic expression of height-diameter relationships. For. Chron. 31, 236–246.

Kienast, F., 1987. Forece: A Forest Succession Model for Southern Central Europe (No. ORNL/TM-10575). Oak Ridge National Lab., TN (USA). Kunstler, G., Coomes, D.A., Canham, C.D., 2009. Size-dependence of growth and mortality influence the shade tolerance of trees in a lowland temperate rain forest. J. Ecol. 97, 685–695. Kunstler, G., Allen, R.B., Coomes, D.A., Canham, C.D., Wright, E.F., 2013. Sustainable management, earthquake disturbances, and transient dynamics: modelling timber harvesting impacts in mixed-species forests. Ann. For. Sci. 70, 287–298. Liu, J., Ashton, P.S., 1995. Individual-based simulation models for forest succession and management. For. Ecol. Manag. 73, 157–175. McGlone, M.S., Hall, G.M.J., Wilmshurst, J.M., 2011. Seasonality in the early Holocene: extending fossil-based estimates with a forest ecosystem process model. Holocene 21, 517–526. Melo, F.P.L., Pinto, S.R.R., Brancalion, P.H.S., Castro, P.S., Rodrigues, R.R., Aronson, J., Tabarelli, M., 2013. Priority setting for scaling-up tropical forest restoration projects: early lessons from the Atlantic Forest Restoration Pact. Environ. Sci. Policy 33, 395–404. Messier, C., Puettmann, K., Chazdon, R., Andersson, K.P., Angers, V.A., Brotons, L., Filotas, E., Tittler, R., Parrott, L., Levin, S.A., 2015. From management to stewardship: viewing forests As complex adaptive systems in an uncertain world. Conserv. Lett. 8, 368–377. Morales, N.S., Perry, G.L.W., Burns, B.R., 2016. Fencing is not enough to reinstate regeneration: evidence from a large fruited canopy tree Beilschmiedia tawa. For. Ecol. Manag. 376, 36–44. Morales, N.S., 2015. Factors affecting recruitment of Beilschmiedia tawa in northern New Zealand. N. Z. J. Bot. 53, 231–240. Ogden, J., Stewart, G.H., 1995. Community dynamics of the New Zealand conifers. In: Enright, N.J., Hill, R.S. (Eds.), In Ecology of the Southern Conifers. Melbourne University Press, Melbourne, pp. 81–119. Perry, G.L.W., Millington, J.D.A., 2008. Spatial modelling of succession-disturbance dynamics in forest ecosystems: concepts and examples. Perspect. Plant Ecol. Evol. Syst. 9, 191–210. Primack, R., Rozzi, R., Feinsinger, P., Dirzo, R., Massardo, F., 2001. Fundamentos de Conservacion Biologica: Perspectivas Latinoamericanas. Fondo de Cultura Economica, Mexico. Reuter, H., Breckling, B., Jopp, F., 2011. Individual-based models. In: Jopp, F., Reuter, H., Breckling, B. (Eds.), Modelling Complex Ecological Dynamics. Springer, Berlin Heidelberg, pp. 163–178. Richardson, S., Smale, M.C., Hurst, J.M., Fitzgerald, N.B., Peltzer, D.A., Allen, R., Bellingham, P.J., McKelvey, P.J., 2009. Large-tree growth and mortality rates in forests of the central North Island, New Zealand. NZES 33, 208–215. Rogers, H.M., 1999. Stand dynamics of Dacrydium cupressinum dominated forest on glacial terraces, south Westland, New Zealand. For. Ecol. Manag. 117, 111–128. Santo-Silva, E.E., Almeida, W.R., Tabarelli, M., Peres, C.A., 2016. Habitat fragmentation and the future structure of tree assemblages in a fragmented Atlantic forest landscape. Plant Ecol. 217, 1129–1140. Scheller, R.M., Mladenoff, D.J., 2007. An ecological classification of forest landscape simulation models: tools and strategies for understanding broad-scale forested ecosystems. Landsc. Ecol. 22, 491–505. Schumacher, S., Bugmann, H., Mladenoff, D.J., 2004. Improving the formulation of tree growth and succession in a spatially explicit landscape model. Ecol. Model. 180, 175–194. Shugart, H.H., Smith, T.M., 1996. A review of forest patch models and their application to global change research. Clim. Change 34, 131–153. Shugart, H.H., 1984. A Theory of Forest Dynamics: the Ecological Implications of Forest Succession Models. Springer-Verlag. Smale, M.C., Bathgate, J.L., Guest, R., 1986. Current prospects for tawa. N. Z J. For. 31, 13–18. Smale, M.C., Dodd, M.B., Burns, B.R., Power, I.L., 2008. Long-term impacts of grazing on indigenous forest remnants on North Island hill country, New Zealand. N. Z. J. Ecol. 32, 57–66. Smale, M.C., Richardson, S.J., Hurst, J.M., 2014. Diameter growth rates of tawa (Beilschmiedia tawa) across the middle North Island, New Zealand–implications for sustainable forest management. N. Z. J. For. Sci. 44, 20. Stewart, G.H., White, J.C., Duncan, R.P., 1998. The structure and dynamics of rimu-dominated forests on glacial moraines, South Westland. N. Z. For. 42, 22–27. Thrippleton, T., Dolos, K., Perry, G.L.W., Groeneveld, J., Reineking, B., 2014. Simulating long-term vegetation dynamics using a forest landscape model: the post-Taupo succession on Mt. Hauhungatahi, North Island, New Zealand. N. Z. J. Ecol. 38, 26–38. West, C.J., 1986. Population Ecology of Beilschmiedia tawa at Pureora Forest. University of Auckland, Auckland, New Zealand, Unpublished PhD Thesis. White, J.W., Rassweiler, A., Samhouri, J.F., Stier, A.C., White, C., 2014. Ecologists should not use statistical significance tests to interpret simulation model results. Oikos 123, 385–388. Wilensky, U., 1999. NetLogo. Center for Connected Learning and Computer-Based Modeling. Northwestern University, Evanston IL.