A large-scale forest landscape model incorporating multi-scale processes and utilizing forest inventory data

Two challenges confronting forest landscape models (FLMs) are how to simulate fine, standscale processes while making large-scale (i.e., .10 ha) simulation possible, and how to take advantage of extensive forest inventory data such as U.S. Forest Inventory and Analysis (FIA) data to initialize and constrain model parameters. We present the LANDIS PRO model that addresses these needs. LANDIS PRO adds density and size mechanisms of resource competition. This is achieved through incorporating number of trees and DBH by species age cohort within each raster cell. Forest change is determined by the interactions of species-, stand-, and landscape-scale processes. Species-scale processes include tree growth, establishment, and mortality. Stand-scale processes include density and size-related resource competition that regulates self-thinning and seedling establishment. Landscape-scale processes include seed dispersal, as well as natural and anthropogenic disturbances. LANDIS PRO is designed to be straightforwardly comparable with forest inventory data, and thus the extensive FIA data can be directly utilized to initialize and constrain model parameters before predicting future forest change. We initialized a large landscape (;10 ha) from historical FIA data (1978) and the predicted forest structure and composition following 30 years of simulation were statistically calibrated against a prior time-series of sequential FIA data (1978 to 2008). The results showed that the initialized conditions realistically represented the historical forest composition and structure at 1978, and the constrained model parameters predicted reasonable outcomes at both landscape and land type scales. The subsequent evaluation of model predictions showed that the predicted forest composition and structure were comparable with old-growth oak forests; predicted forest successional trajectories were consistent with the expected successional patterns in oak-dominated forests in the study region; and the predicted stand development patterns were in agreement with the established theories of forest stand development. This study demonstrated a framework for forest landscape modeling including model initialization, calibration, and evaluation of predictions.


INTRODUCTION
Forest changes as the result of stand-scale processes such as resource competition, forest landscape processes (FLPs) such as fire and harvest, and regional processes such as climate change.Forest landscape models (FLMs) have increasingly become important tools for predicting these changes over large spatial scales (landscape and regional) (Keane et al. 2004, Perry and Enright 2006, He 2008, Taylor et al. 2009).There has been remarkable progress over the past two decades in the development of FLMs.FLMs have been used to examine the effects of wildfire (e.g., Yang et al. 2004, Keane et al. 2008), prescribed fire (e.g., Syphard et al. 2011), harvest (e.g., Gustafson et al. 2000), insect outbreaks (e.g., Sturtevant et al. 2004), and climate change (e.g., Schumacher and Bugmann 2006, Gustafson et al. 2010, Thompson et al. 2011) on forest changes at landscape scales.Despite their usefulness, a fundamental challenge remains: how to simulate fine, stand-scale processes while making large-scale simulation possible.
FLMs usually treat the landscape as raster cells.Site-scale processes such as succession and competition are simulated in each cell and FLPs are simulated over a subset of spatially continuous cells.Site-scale processes include speciesand stand-scale processes.Species-scale processes include tree growth, establishment, and mortality.Stand-scale processes include competition for resources (e.g., light, water, and nutrients) among trees.Stand-scale processes regulate species-scale processes, such as competition-caused mortality (Oliver and Larson 1996).FLPs also interact with species-scale processes through disturbance-induced mortality.Because the computation loads for FLMs increase exponentially with increasing site-scale complexity and size of landscape simulated, FLMs typically simplify species-and stand-scale processes in order to simulate a relatively large landscape (Mladenoff 2004).Tree information in early FLMs at site-scales is aggregated as absence/presence of species age cohorts in LANDIS (Mladenoff and He 1999) or forest type in LANDSUM (Keane et al. 2002) and SIMPPLLE (Chew et al.2004) in lieu of tree density and size.Tree-growth in many FLMs is simplified as age increment (e.g., Li and Barclay 2001) rather than DBH and crown increment as a function of resource competition, particularly light in gap models (Pacala et al. 1996, Deutschman et al. 1999).Resource competition at stand scales is either simplified using a species vital attributes approach in which shadetolerant species outcompete shade intolerant species (e.g., LANDIS), or using a succession pathways approach in which late-successional forest types replace early-successional forest types through predefined transition probabilities (Chew et al. 2004, Keane et al. 2004).
The simplifications of site-scale processes in the early FLMs make them applicable to landscapes where site-scale processes can be strongly influenced or overridden by FLPs (e.g., clearcutting or stand replacing fire) (Romme et al. 1998, Schumacher et al. 2004, Turner 2010).However, the simplified FLMs lack primary mechanisms of stand-scale processes such as density and size related competition (Moorcroft et al. 2001, Purves et al. 2008).Thus, they may not be adequate for landscapes where the effects of disturbance are relatively weak and site-scale processes are critical determinants of future forest changes (Schumacher et al. 2004, Bohlman andPacala 2012).
Recent FLMs have advanced in incorporating species-scale quantitative information (e.g., tree number, biomass), thus improving the mechanisms of stand-scale resource competition.For example, in LANDIS II, biomass for each species age cohort is added and a ratio of actual biomass to potential biomass for a cell is used to quantify resource availability (growing space), assuming species-age cohort biomass implicitly incorporates density information (Scheller et al. 2007).Potential biomass in LANDIS II is derived empirically (Smith et al. 2006) or assumed as 30-fold of maximum aboveground net primary productivity (ANPP) for each species, which is estimated from an ecophysiography model PnET II (Aber et al. 1995).Besides simulating biomass dynamics, LANDIS II implements ecosystem processes such as nutrient cycling in each cell, which expands the scope of traditional FLMs (Scheller et al. 2011).LANDCLIM tracks number of trees and biomass by species age cohort.It is unique in that it introduces gap model dynamics into simulating site-scale dynamics by including the interactions of abiotic variables (terrain, soil, and climate) and biotic variables (tree size, density, and biomass).Stand-scale resource competition is simulated by growth-and densitydependent mortality determined by maximum stand biomass (Schumacher et al. 2004).Because of the additional computation loads required for simulating the greater complexities at each cell, the reported simulation capacity of these models is moderate, for example, 10 4 -10 5 cells in Schumacher and Bugmann (2006); 10 4 cells in Scheller et al. (2011), 10 6 cells in Thompson et al. (2011).Thus, in terms of simulation capacity, these models have not demonstrated substantial improvement compared to a decade ago, for example, 10 6 cells in He and Mladenoff (1999a).
Another challenge confronting FLMs is that they have not fully taken advantage of extensive forest inventory data to initialize and constrain model parameters.This is because model input parameters such as absence/ presence, biomass, and ANPP are not straightforwardly comparable with forest inventory data that usually consists of the number of trees and stem diameters (Pacala et al. 1996, Moorcroft et al. 2001, He et al. 2011).Currently, most FLMs only utilize tree species occurrence and frequency imputed across the landscape (e.g., Zhang et al. 2009, Thompson et al. 2011).Meanwhile, millions of records of individual trees are now available in the U.S. Forest Service Forest Inventory and Analysis database (Woodall et al. 2010) and elsewhere, which provide a new opportunity for FLMs to directly apply these data to initialize and constrain model parameters to improve the future model predictions using data assimilation (DA) techniques (Luo et al. 2011, Peng et al. 2011).DA techniques integrate inventory data with ecological models to constrain the initial conditions and parameters; thus, the simulated results can best match the observed before applying models to future predictions (Williams et al. 2005).
All of the above add up to the need for innovative approaches for FLMs to incorporate key mechanisms that determine stand-scale resource competition while making large-scale simulation feasible, and to directly employ forest inventory data for model initialization and constraining parameters.Individual-based models (IBMs) address many of the issues discussed above and are comparatively easy to parameterize and calibrate against inventory data (Pacala et al. 1996).Applying IBMs to spatially explicit simulation over a large scale, however, would require immense computation loads (Moorcroft et al. 2001, Urban 2005, Seidl et al. 2012).To solve this problem, what is needed is to scale up individuals to retain only individual-level mechanisms without simulating every tree.Recent progress in IBMs to scale up individuals by using mean density and size structure of species offers a possible resolution to tackle the challenges confronting FLMs (Strigul et al. 2008).Once the density and size information is available, well tested theories in forest dynamics such as stand density index (SDI) (Reineke 1933) and selfthinning (Yoda et al. 1963) can be readily applied to quantifying stand-scale resource competition and simulating competition-caused mortality.
We present an innovative approach implemented in LANDIS PRO, a variant from LANDIS FLM family.We demonstrate how forest change is simulated by incorporating species-, stand-, and landscape-scale processes; individual tree information is scaled up through tracking the number of trees by species age cohort.Standscale resource competition is quantified by the amount of growing space occupied (GSO).Species-scale processes including growth, establishment, and mortality are regulated by standscale resource competition as a function of environmental heterogeneity and resource availability.These processes, in combination with landscape disturbance and seed dispersal largely determine the dynamics of forest composition and structure.We show how LANDIS PRO initialization and predictions are straightforwardly comparable with forest inventory data, and thus, the model can be directly initialized from FIA data, as well as calibrated against FIA data.Specifically, we demonstrate that the model can be directly initialized from a historical FIA inventory (1978), calibrated against a prior timeseries of sequential FIA inventories (1978 to 2008).We also demonstrate that the model can be applied on a forest landscape 10 7 ha in spatial extent.Finally, we show how model predictions are evaluated using stand density management diagrams (Gingrich stocking charts and Reineke density diagrams) and the existing empirical studies.

METHODS
The overall design and structure of LANDIS PRO LANDIS PRO is a raster-based FLM that evolved over 15 years of development and applications of the original LANDIS model (Mladenoff et al. 1996, He and Mladenoff 1999a, Mladenoff 2004, Yang et al. 2011).LANDIS PRO simulates forest change over large spatial (;10 8 ha) and temporal (;10 3 years) extents with flexible spatial (10-500 m pixel size) and temporal resolutions (1-10 years).In LANDIS PRO, forest succession and dynamics are designed to incorporate species-, stand-, and landscape-scale processes (Fig. 1).Species-and stand-scale processes are simulated within each raster cell and landscape-scale processes are simulated across the whole landscape, and these three scales processes interact with each other.Within each raster cell, the model records number of trees by species age cohort (Fig. 1) and size (e.g., mean DBH) for each age cohort is derived from empirical age-DBH relationships (e.g., Loewenstein et al. 2000).Thus, number of trees by species age cohort is the only parameter for the landscape initialization.With available density and size information, LANDIS PRO key stand parameters including density, basal area, stocking, biomass, carbon, and importance value for each species at each raster cell.Taken together, LANDIS PRO initialization and simulation results allow for straightforward comparisons with forest inventory data, thus making it possible for LANDIS PRO to fully utilize current extensive forest inventory data to constrain initial forest conditions and model parameters.
LANDIS PRO stratifies the heterogeneous landscape into land types (also called ecoregions for broad-scale studies), which capture the coarsest-level (coarse grain) spatial heterogeneity in resource availabilities (maximum growing space occupied) and species assemblages (species establishment probability (SEP)).SEP is a number from 0.0-1.0,reflecting how different environmental conditions favor a particular species in terms of its seedling establishment.SEP is the same parameter used in previous LANDIS models (Mladenoff and He 1999), LANDIS II (Scheller et al. 2007), andLANDCLIM (Schumacher et al. 2004), which are estimated from a coupled landscape and ecosystem modeling approach to simulate the effects of climate change on seedling establishment.For example, SEPs used in LANDIS PRO are derived from LINKAGES II (Wullschleger et al. 2003)

Landscape-scale process: seed dispersal
Landscape-scale processes simulated in LANDIS PRO include seed dispersal (exotic species invasion), fire, wind, insect, disease spread, forest harvesting, fuel treatments, and silvicultural treatments.These disturbance processes release growing space on one or more stands on the landscape, thus resetting the development stage of affected stands (Oliver and Larson 1996).In this paper, we focus only on species-, stand-scale processes, and landscapescale seed dispersal.The simulation of other landscape-scale processes such as harvest and fire in LANDIS PRO is described elsewhere (e.g., Fraser et al., in press).
Seed dispersal is a spatial process in LANDIS PRO.It is simulated using a dispersal kernel determined by species-specific maximum dispersal distances, where the probability of seed dispersal to every cell is calculated using a negative exponential decay function (He and Mladenoff 1999b).The number of potential germination seeds (NPGS) reaching a cell from the parent tree is the product of dispersal probability and number of potential germination seeds.The total NPGS for each species reaching a given cell is accumulated from all available mature trees within the dispersal kernel.NPGS for each species is a new parameter added to LANDIS PRO that influences species density and basal area.There is limited information regarding this parameter.Burns and Honkala (1990) provided some information about this parameter on which further calibration is based to ensure predicted density and basal area match with forest inventory data.Thus, LANDIS PRO accounts for dispersal limitation and seed availability that constrain the potential species distriv www.esajournals.orgbution under rapidly changing environment.

Stand-scale processes
Stand-scale processes include resources competition that regulates competition-caused mortality (self-thinning) and seedling establishment.The competition intensity is quantified by the amount of growing space occupied (GSO) within each cell.The self-thinning is implemented using Yoda's self-thinning (Yoda et al. 1963), and seedling establishment is determined by available growing space, species shade tolerance, and SEP.Resource availability varies among different stand development stages because of the dynamics of establishment and mortality.The simulation of stand development patterns is governed by GSO and follows well documented stages of stand development.
Resource competition quantified by growing space occupied (GSO).-TheReineke stand density index (SDI) is used to quantify stand-scale resource competition in pure stands by using the quadratic mean stand diameter and number of trees (Reineke 1933, Long and Daniel 1990, Shaw 2000).However, estimation of SDI in mixed stands has been an ongoing research because of the differences in growing space requirements among species and variation caused by uneven or irregular stand age structures (Woodall et al. 2003).We use GSO to quantify stand-scale resource competition in LANDIS PRO to avoid accounting for the differences in growing space required by different species.We estimate GSO by summing the total minimum growing space required to support to all trees on a site.The minimum growing space required for each tree at a given species and size is derived from the Reineke stand density index (SDI) and maximum SDI.Maximum SDI is defined as the maximum number of 10-inch diameter trees per hectare and has already been reported for many species (Reineke 1933, Long 1985).
We calculate GSO by converting the number of trees by species age cohort (DBH) (Fig. 2a) into the number of 10-inch diameter trees, which we call number of standard trees (NST i ), by using Reineke SDI (Fig. 2b, Eq. 1); We calculate v www.esajournals.orgminimum growing space required for each 10inch diameter tree by species (MinGS i _standard) using the inverse of Maximum SDI for given species (Eq.2; Reineke 1933); the total growing space is the sum of total minimum growing space occupied by all 10-inch diameter trees.The GSO for the cell is calculated as the percentage of the cell area occupied by all 10-inch diameter trees (Fig. 2c,Eq. 3).Since this percentage is a summation of values associated only with tree DBH and density, GSO is independent of site quality and can be compared for stands with mixed species and age classes as well as for evenaged monocultures.
where DBHi ij (in cm) and NT ij (in stems) are the mean diameter and number of trees for jth age cohort of species i; MaxSDI i is the maximum number of trees with 10-inch (25.4-cm) diameter per hectare for species i; raster_area represents the area for each raster in the model in m 2 .Stand development patterns.-Standdevelopment patterns follows (1) stand initiation stage, (2) stem exclusion stage, (3) understory reinitiation stage, and (4) old-growth stage (Oliver and Larson 1996).Four specific thresholds of GSO are user-defined by land type to regulate seedling establishment in stand initiation stage, where seedlings can only become established before stands reach fully occupied (Fig. 3, Table 1): (1) open grown (0-GSO 1 ), (2) partially occupied (GSO 1 -GSO 2 ), (3) crown closure (GSO 2 -GSO 3 ), (4) fully occupied (GSO 3 -MGSO).MGSO represents the maximum growing space that can be occupied.Once stands exceed MGSO, stands are presumed to reach stem exclusion stage.Meanwhile, self-thinning is initiated and continues to the following understory re-initiation and old-growth stages (Oliver and Larson 1996).
An open, recently disturbed forest stand enters the stand initiation stage, which is characterized by widespread establishment of primarily shade intolerant species and rapid growth of advance reproduction (Franklin et al. 2002).As the stand is progressively filled and occupied, only seedlings of shade-tolerant species can become established.When the growing space becomes fully occupied through a combination of tree establishment and growth, the stand is considered fully stocked, thereby, entering the stem exclusion stage of development.Trees that are small, shade intolerant, or approaching their species' longevity can be outcompeted first via self-thinning (Reynolds andFord 2005, Coomes andAllen 2007).As the mean size of trees in the stand increases, larger canopy gaps are created by tree mortality.During the understory reinitiation stage, these gaps are refilled by establishment of new seedlings or the lateral growth of adjacent of trees.Continued tree growth and mortality in the absence of exogenous disturbance moves the stand into the oldgrowth stage of development.At old-growth stage, old trees die as they reach their species' longevity, creating large canopy gaps that promote tree regeneration and move the stand into an uneven-aged condition.
The definitions of the four GSO thresholds can vary by landtype and can be modeled outside LANDIS PRO platform to incorporate the effects of varying MGSO resulting from nitrogen or CO 2 fertilization under climate change.They can also be defined for a variety of ecosystems.For example, a savanna system may only have first development stage and have a low GSO 1 threshold; and a woodland system may never reach the crown closure stage and may have low GSO 1 and GSO 2; MGSO thresholds for southfacing land types may be lower than those for north-facing land types because of the higher moisture and nutrient availability (Kabrick et al. 2008, Johnson et al. 2009).
Self-thinning.-We simulate resource competition-caused tree mortality in LANDIS PRO as self-thinning, which is a natural process because of limited growing space and less shade tolerance (Monserud et al. 2005, Shaw 2006).The selfthinning associated tree mortality is characterized by a decrease in the number of trees with increasing average tree size in the stand and follows the À3/2 rule (also referred as Yoda's selfv www.esajournals.orgthinning line) (Yoda et al. 1963).
In LANDIS PRO, MGSO is defined to identify the trajectory of the self-thinning line.When stands approach or exceed MGSO (the upper limit of the available growing space), the model initiates stand-scale self-thinning.Specifically, four quarterly age ranges for thinning are created based on species longevity.Self-thinning will be performed progressively through these four age ranges until sites converge with self-thinning line.In each age range, growing space occupied for one standard tree by species age cohort first ranks based on the growing space occupied.Following the rank, the mortality by species age cohort is determined as a function of age and shade tolerance of the species.In general, higher shade tolerance species have lower thinning percentage, because they are more shade tolerant and can survive longer period of suppression (Reynolds and Ford 2005).Thus, stand development trajectories converge with the self-thinning line and move along the self-thing line from lower right to upper left (Fig. 3).

Species-scale processes
Species-scale processes in LANDIS PRO include growth, seedling establishment, sprouting, and sexual mortality.We simulate species-scale processes using species vital attributes including species age of maturity, longevity, shade tolerance, maximum DBH, and average number of seeds per mature tree per year, and Maximum SDI (Table 2).We simulate tree growth as age and DBH increment by species age cohort.Age increment is determined by the model time step.DBH increment is simulated using empirical log v www.esajournals.orgnormal distribution (Condit et al. 2006) or calibrated locally by land type (Murphy 1998, Loewenstein et al. 2000), which can account for the effects of varying resource availability by land type and climate change on tree growth.
Seedlings in LANDIS PRO include seed variability and sprouting.Sprouting is simulated for trees with capability for vegetative reproduction following exogenous disturbances such as harvesting, wind, and fire.Species-specific sprouting probabilities are age dependent and specified as part of each species' vital attribute data (Table 2).Seedling establishment at stand initiation stage is regulated by the NPGS, GSO, species shade tolerance, and SEP (Table 1).
Tree mortality in the model includes: (1) natural mortality due to reaching species longevity, (2) background mortality caused by unknown or unspecified factors (Hamilton and Edwards 1976), (3) stand-scale resource competition-caused mortality (self-thinning), and (4) landscape disturbance and management-caused mortality.Compared to the previous versions of LANDIS, the latter two forms of mortality simulation are new and are tree-based as opposed to age cohort-based (i.e., rather than modeling the death of an entire age cohort as a single entity, individual trees or subsets of trees within a cohort can die).Tree-based mortality improves the simulation realism for disturbance and competition, which often remove only part of age cohort.The background mortality allows for customizing an empirical function estimating the mortality caused by exogenous factors such as climate change (Hamilton andEdwards 1976, Johnson et al. 2009).

Study area
The study area was over 10 7 ha and comprised of 4100 3 3200 cells with a resolution of 90 m (0.81 ha).Boundaries corresponded to FIA Survey Unit 5 in Arkansas (Fig. 4).This area encompasses the Ozark and Boston Mountains in Arkansas.It is characterized as deeply dissected and rugged, with elevations ranging from 213 m in valley bottoms to 762 m at the highest ridge crests.The average annual temperature and precipitation range from 148 to 178C, and from 1150 to 1325 mm, respectively; the majority of precipitation occurs in the spring and fall.The area is characterized as mixed hardwood-pine forest dominated by oak (Quercus spp.), hickory (Carya spp.), and shortleaf pine (Pinus echinata Mill.).The dominant hardwood species include white oak (Q.alba L.), post oak (Q.stellata Wangenh.), chinkapin oak (Q.muehlenbergii Engelm.),black oak (Q.velutina Lam.), northern red oak (Q.rubra L.), blackjack oak (Q.marilandica Muenchh.), southern red oak (Q.falcate Michx.), pignut hickory (C.glabra Sweet), and black hickory (C.texana Buckl.).Species composition and distribution were significantly altered since European settlement (Chapman et al. 2006).A large portion of the current forest regenerated after extensive timber harvest in the early 1900s and today the age of dominant and codominant oaks typically ranges from 60 to 90 years.The stem density has greatly increased reaching full stocking because of nearly a century of fire suppression (Heitzman 2003).
Because prior disturbances were minimal in this area (less than 5% of FIA plots recorded disturbance) over the last 30 years, we only used undisturbed FIA plots to initialize the landscape and we only simulated forest growth and  and (2) no evidence of disturbance including logging, insects, disease and fire since the last measurement.We estimated tree ages for each FIA plot using an empirical age-DBH equation (Loewenstein et al. 2000) and aggregated the number of trees and basal area by species age cohort per plot and per hectare based on tree expansion factors recorded in FIA (Woodall et al. 2011).Stratum expansion factor, an area expansion factor in FIA was used to scale-up FIA plot inventory to the land type or landscape scale to estimate number of trees and basal area for each species by land type and the whole landscape.

Landscape initialization
Species vital attributes and land types.-Wederived species vital attributes (Table 2) from existing data sets for the Boston Mountains (Spetich and He 2008) and Silvics of North America (Burns and Honkala 1990).We grouped tree species into six functional groups accounting for over 90% of total basal area: white oak (white oak and post oak), red oak (northern red oak and southern red oak), black oak, hickory, pine (shortleaf pine and loblolly pine, Pinus taeda L.), and maple (red maple, Acer rubrum L. and sugar maple, A. saccharum Marsh.).
The land type map including five land types was derived based on slope position and aspect from the Ozark-St.Francis National Forest Ecological Classification System (ECS) and Landtype Associations (LTAs) (Fig. 4).SEPs by species for each land type were derived from existing forest inventory data for the Boston Mountains, Arkansas (Spetich and He 2008) that were estimated from ecosystem model LINKAGES II (Wullschleger et al. 2003).Four values of GSO (GSO 1 , GSO 2 , GSO 3 and MGSO) were defined by land type and calibrated against FIA data.
The initial forest conditions.-Theinitial forest species composition map (1978) for the study area containing number of trees by species age cohort in each cell was created from the 1978 FIA data (Fig. 5a) using Landscape Builder (Dijak, in press).This program first converted the stemlevel inventory data within each FIA plot to number of trees per hectare by species age cohort using published DBH-age relationships.These plot-level data were then scaled for individual cells based on the cell size and FIA tree expansion factor (Woodall et al. 2011).Finally, Landscape Builder stochastically selected and assigned a representative FIA plot to each cell to represent species composition in the initial landscape.Each FIA plot assigned to represent a particular cell was screened and stratified to draw only from the pool of FIA plots that matched the cell in FIA unit, national forest type, and national forest size class.

Constraining the initialized forest conditions and model parameters
To constrain the initialized landscape to ensure it represented the historical forest conditions, we iteratively adjusted species growth rates until the initialized species basal area matched with the summarized FIA data of 1978 at both landscape and land type scales (pass chi-square test for no differences in species density and basal area).We then used the initial landscape for 1978 as the starting point and simulated forest succession to 2008 (30 years) to calibrate the NPGS for each species.This was done by iteratively comparing the simulated number of trees and basal area by v www.esajournals.orgspecies against the observed changes in the FIA inventory for same time period at landscape and land type scales (pass chi-square test for no differences in species density and basal area by land type).

Evaluation of model predictions
We used the constrained model to simulate forest change from 1978 to 2128 (150 years).Ideally evaluation of model predictions requires independent data (Clark et al. 2001, Gardner and Urban 2003, Moorcroft 2006), however, independent data for comparing with model predictions are not often available for a landscape scale study (Shifley et al. 2009).Thus, we compared our model predictions with field studies of upland Missouri old-growth forests (Richards et al. 1995, Shifley et al. 1995) and with the established theories of forest succession and stand development.Four priors in comparing with field studies were (1) whether high mortality was observed in the early simulation stage, (2) whether the predicted maximum basal areas in the later simulation stage was consistent with the basal areas measured for mature, old-growth oak forest in this region (Shifley et al. 1995, Richards et al. 1995), (3) whether the oak-dominated forests will be gradually successionally replaced by longer-lived white oak species and shadetolerant species such as maple in absence of disturbance, (Johnson et al. 2009), and (4) whether the predicted total density, basal area, biomass, and carbon for northeast land type were higher than for southwest land types, because, in our study, southwest and northeast land types presented large differences in physical environments.
We used Gingrich (1967) stocking charts and Reineke (1933) density diagrams to evaluate stand development patterns, specifically the relationship between basal area, density, and quadratic mean diameter over time.Gingrich (1967) stocking charts for upland oak forests indicated combinations of mean stand diameter, basal area per acre, and number of trees per acre at which a forest stand was considered fully stocked (i.e., capable of fully occupying but not exceeding the available growing space) (Ernst and Knapp 1985).The upper limit of stand occupancy was indicated by the line for 100% stocking (often termed the A-line).The minimum conditions at which the trees on the site can fully occupy the growing space occur at approximately 58% stocking (often termed the B-line).In theory, an undisturbed oak-hickory stand at a stocking level less than 100% would gradually increase in basal area and decrease in number of trees at rates that would move the stand toward but not consistently above 100% stocking.Shifley et al. (1995) demonstrated that the stocking percentage for mature and old-growth upland oak forests in Ozark Highlands was within 80% ; 95%.Reineke (1933) density diagrams, which were algebraically analogous to the Gingrich stocking guides, provided another graphical framework to examine trajectories of mean stand conditions over time with respect to available growing space.

Constraining the initial forest conditions
Forest composition for the initial forest conditions.-Theinitial landscape constructed with Landscape Builder from FIA data for 1978 captured the historical species composition of oak-dominated forests at 1978 reasonably well.There was no significant difference in species density (v 2 ¼ 1.93, df ¼ 5, P ¼ 0.86) or basal area (v 2 ¼ 1.40, df ¼ 5, P ¼ 0.92) at the landscape scale nor by land type (southwest: 6).The white oak group (consisting of white oak and post oak) was the predominant species group across the landscape, comprising 35% of the total basal area.The red oak group (northern red oak, southern red oak) and black oak together included another 30% of the total basal area.Hickory was consistently abundant across the landscape, making up 20% of the basal area.Pine and maple followed in abundance with 10% and 5% of the total basal area, respectively.
Forest structure for the initial forest conditions.-Therewas no significant difference between the initialized 1978 forest structure and the FIA data at the landscape scale for density by age classes (v 2 ¼ 3.23, df ¼ 9, P ¼ 0.95) or basal area by age classes (v 2 ¼ 1.07, df ¼ 9, P ¼ 0.98) for all species combined (Fig. 7).This illustrated that the initial forest landscape was capable of capturing much of the variation in forest structure.The age class distribution (estimated from the DBH distribution) for each species had a reverse J-shape both for the initialized landscape and for the FIA data.In 1978, species in the white oak and the red oak groups were dominant in the overstory, hickory was of intermediate abundance, and maple was predominantly an understory species group.

Evaluation of model predictions
Forest composition and structure.-Asexpected, the density of white oak, red oak, black oak, hickory and pine decreased dramatically over the 150-year period as a consequence of self-thinning associated with forest maturation and tree growth.Simulated basal area, biomass and carbon increased from 1978 to a peak at 2098, followed subsequently by slight declines, because a large proportion of trees in the red oak and black oak species groups that became established in the early to mid1900s reached longevity, died, and were replaced by young trees.
Predicted basal area reached a maximum of 100 ft 2 /acre on southwest land types and 120 ft 2 / acre on northeast land types.These predicted values were consistent with the basal area estimates of 102-120 ft 2 /acre reported by Shifley et al. (1995) and Richards et al. (1995) for mature, undisturbed forests in the Ozark Highlands.Also, in year 2128 (Fig. 5c), species basal area, biomass, and carbon were concentrated in large, old trees; this pattern was associated with the shift in age distribution due to continued diameter growth of trees in the absence of exogenous disturbances.In year 2128 the age distribution of trees greater than 96 years old approaches a bell-shape with ten or fewer trees per acre in younger diameter classes (Fig. 7).In upland oak forests in the Central Hardwood region, competitive oak regeneration is highly dependent on disturbance (e.g., fire and forest harvesting).Without disturbance, oak-dominated forests frequently are successionally replaced by a gradually increasing population of shadetolerant species (e.g., sugar maple, sweet gum) (Fig. 8).This is one of the most pervasive problems associated with sustaining oak-dominated forests, particularly on mesic sites (Spetich 2004, Dey et al. 2010).
Successional dynamics.-Ourmodel predictions indicated that without disturbances white oaks would continue to dominate the landscape for v www.esajournals.orgthe next 150 years (Fig. 5c, Fig. 8).The red oaks and especially black oak declined in basal area after 2080 as many trees of those species groups reached their maximum longevity and died.The lack of simulated disturbance favored establishment of shade-tolerant maple species, thus the predicted maple density gradually increased from 5% in 1978 to 20% in 2128.These predicted forest successional trajectories were consistent with our expectations for oak-dominated forests in the study region: in absence of disturbance, mixed oak forests typically transition to a greater proportion of longer-lived white oak species and shade-tolerant species such as maple increase in abundance (Johnson et al. 2009).
We plotted the stand development trajectories from 1978 to 2128 for mean stand conditions for the entire landscape and by land type on Gingrich stocking charts (Fig. 9a).Trajectories remained within the fully stocked zone (between 58 and 100%).Mean stocking percentages increased over time, as would be expected without exogenous disturbances.Also as expected, mean stocking percentages on northeast land types increased more rapidly and the mean number of trees decreased more rapidly than on southwest land types that typically had lower site quality and slower tree growth.At year 2128 of the simulation period, the mean stocking percentage reached 90% for northeast land types and 75% for southwest land types.These predicted stocking percentages were consistent with the stocking percentages of 80%-95% reported by Shifley et al. (1995) for mature and old-growth upland oak forests in the adjacent Ozark Highlands.We also illustrated the increase in mean diameter and decrease in number of trees per acre associated with self-thinning as stands grow undisturbed at landscape and land type scales with the Reineke density diagrams (Fig. 9b).
Landscape heterogeneity.-We were able to capture landscape heterogeneity among land types with our modeling approach.Southwest land types and northeast land types presented large differences in physical environments which were reflected in differential rates of forest change among land types.Our predicted results showed that the total density, basal area, biomass, and carbon for northeast land types were on average higher than those for southwest land types (Fig. 8).The maximum stocking percentage reached at northeast land type was also higher than for southwest land types (Fig. 9a).Therefore, user-defined parameters were able to capture the effects of the landscape heterogeneity in species composition, biomass, and carbon.

Design implications
The approach of incorporating stand-scale processes in FLMs has benefited from the recent advances in IBMs that scale up individuals through simulating number of trees and size by species age cohort without simulating every tree (Moorcroft et al. 2001, Strigul et al. 2008).Stand- v www.esajournals.orgscale resource competition is quantified by GSO using density and size information.Such information is the primary determinant of resource competition deriving forest dynamics (Moorcroft et al. 2001, Purves andPacala 2008).Resource competition regulates seedling establishment and self-thinning, accounting for resource availability at different stand development stages.Our approach in LANDIS PRO directly simulates the competition-caused mortality using À3/2 power rule of self-thinning (Yoda et al. 1963), which explicitly characterizes the decrease in the number of trees with increasing average tree size in the thinning stand.
Our intent was to find a balance between stand-scale complexity and computation capacity.We demonstrated that the LANDIS PRO makes large-scale simulation (;10 7 ha) possible while adding stand-scale complexity.Computation efficiency in our approach was achieved by implementing rather than simulating the emergent properties of stand development, analogous to the concepts of microscopic equations used to scale up individual tree level information to stand scale (Strigul et al. 2008).Simulating forest change at large scales is necessarily to incorporate great environmental heterogeneity.For example, Yang et al. (2011) demonstrated that species abundance and composition differed significantly when changing the spatial extent of simulation from 10 3 3 10 3 cells to 10 4 3 10 4 cells with 30 m resolution in Missouri central hardwood forests.Simulating at large spatial scale is also valuable for comparing FLM predictions to those made from niche and (biogeochemical and biogeographical) process models.Currently, niche and process models are the primary tools for predicting the future climate effects on forest composition and distribution at regional scales (Beaumont et al. 2007, Morin and Thuiller 2009, Scheiter and Higgins 2008), whereas FLMs have contributed little to this type of research due to their incapability to simulate sufficiently large landscapes.Recent studies have shown that FLPs v www.esajournals.orgsuch as timber harvest may exert greater effects on forest change than direct change caused by climate warming (Schumacher andBugmann 2006, Gustafson et al. 2010).The improved simulation capacity may allow for incorporating the effects of disturbance and management at regional scales thus reducing prediction uncertainties (Purves and Pacala 2008, Iverson et al. 2011, Matthews et al. 2011).

Data assimilation and the evaluation of model predictions
Our incorporation of density and size makes LANDIS PRO easy to parameterize and compare with field data, which enables us to incorporate DA techniques to constrain the initial conditions and model parameters against FIA data.Thus, the predicted results match the observed data as closely as possible before predicting the future changes.The DA procedures are increasingly important in improving ecological forecasting and estimates of uncertainty (Williams et al. 2005, Hobbs and Ogle 2011, Luo et al. 2011).Advanced capability for predicting future changes is considerably necessary, because the success in mitigating and adapting to future changes in disturbance and climate rests with our capacity to predict (Clark et al. 2001, Moorcroft 2006, Clark and Gelfand 2006, Coreau et al. 2009).
We have shown that the initialized 1978 forest conditions from FIA data reasonably represented the historic forest structure and composition and that predicted forest structure and composition following 30 years of simulation could be statistically calibrated against FIA data.The FIA data have a relatively short time span (e.g., 1978-2008 for our study area), and thus temporal autocorrelation (Arau ´jo et al. 2005) may limit the effective use of FIA data to validate predicted vegetation change for long projection periods.Nevertheless, FIA data provide rare spatial time series that cover a wide area from which various forest successional stages across space may mediate the relatively short survey time span.Finding truly independent large-scale calibration and validation data sets is problematic.However, over time longer series of repeated FIA measurements will accumulate and at least partially alleviate this concern.Likewise, a limited number of landscape-scale experiments that measure tree-, stand-, and landscape-scale changes over time (e.g., Shifley andKabrick 2002, Hardwood Ecosystem Experiment 2012), will gradually provide suitable independent data for calibrating and evaluating landscape model predictions over multiple decades towards a model-data infusion approach, an emerging area of research in ecology (Moorcroft 2006, Clark and Gelfand 2006, Coreau et al. 2009, Luo et al. 2011, Peng et al. 2011).
Evaluating the model predictions is a critical process in ecological modeling to quantify the reliability and confidence of model predictions (Rykiel 1996, Clark et al. 2001, Shifley et al. 2009).The evaluation results showed that the predicted succession trajectories followed the patterns generally observed in central hardwood forests (Abrams 2003, Johnson et al. 2009, Dey et al. 2010): with the absence of disturbance, oakdominated forests were predicted to be successionally replaced by mesic, shade-tolerant species (e.g., maples).Even though the model is more mechanistic compared to previous versions of LANDIS, it uses empirical relationships to simulate tree growth, mortality, and stand development.Such empirical relationships need verification.Empirical studies of old-growth oak forests have been used to evaluate whether the predicted results are ecologically reasonable.Gingrich stocking charts and Reineke density diagrams have been used to evaluate the model design including stand development and patterns of succession.The overall model evaluation showed that the model realistically predicts patterns of succession and old-growth forest structure and composition.Simultaneously, the underlying model theories and design have been evaluated.Through this study we have demonstrated a more thorough framework for FLM simulation including model initialization, constraining the initial landscape and model parameter that have been relatively less explored for FLMs applied to over such a large spatial extent.

Model applications and future research
LANDIS PRO has several potential applications.First, the model provides quantitative information (e.g., basal area, density, biomass, carbon, important value), which can be applied to address the primarily critical questions regarding forest composition, structure, biomass, carbon, and species biodiversity.These questions v www.esajournals.orgincreasingly concern scientists and society under the future climate change and ecosystem management (Thuiller 2007, Frelich andReich 2009).Second, by explicitly incorporating the effects of climatic factor on seedling establishment and background mortality, the model can be suitable for addressing the impacts of climate change on forest dynamics.Third, by simulating the number of trees and DBH by species age cohort, LANDIS PRO can improve in simulating natural disturbance and forest management scenarios, which greatly affect forest structure and dynamics.For example, varying fire severity can be simulated by modeling the proportion of killed trees by species age cohort.Harvest can be simulated in greater details with silvicultural prescriptions based on density, basal area, or stocking percent of the current and/or the residual stand (Fraser et al., in press).Fourth, because of the greater level of species-and standscale details carried in the model, standard standscale measurements (e.g., Reineke density diagram and Gingrich stocking chart) can be used to analyze the simulation results.This makes model predictions more relevant to forest management and planning.
Our additions to LANDIS PRO of competitioncaused mortality, of stocking control per unit area, and of growing space-regulated seedling establishment continue the evolution of FLMs that incorporate greater detail and realism at stand scales.Limits on computation capacity and data for model development and landscape initialization have constrained faster progresses in integration of species-, stand-, and landscapescale processes in one simulation model.Nevertheless, we envision the day when a future variant of FLMs could model stand-scale dynamics in much greater detail.For example a future LANDIS model might simulate tree and standscale forest change by invoking an established individual-tree-based simulation model such as the Forest Vegetation Simulator (Dixon 2002) for each raster at each time step.Forest growth and yield modelers who work primarily at the tree and stand scale have been simultaneously striving to incorporate landscape-scale processes within their modeling systems (e.g., Crookston and Stage 1991, Falkowski et al. 2010, Landscape Management System 2012).Further integration of existing tree-, stand-, and landscape-scale modeling expertise seems inevitable and highly desirable for addressing complex forest management issues.Although landscape-scale forest modelers have often worked independently from tree-or stand-scale modelers (e.g., those focused on timber growth and yield; see Weiskittel et al. 2011), greater collaboration is likely to be beneficial in developing robust, multi-scale models.

Fig. 2 .
Fig. 2. Procedures used to estimate growing space occupied (GSO) using Reineke stand density index (SDI) and maximum SDI in LANDIS PRO.

Fig. 4 .
Fig. 4. The 10 7 ha study area is located in Northern Arkansas within FIA survey unit 5 as indicated by the green area on the top panel.The study area is dominated by oak forest, with a variety of land types.

Fig. 6 .
Fig. 6.Comparison of simulated density and basal area by species group against FIA data at 1978 and 2008 at landscape and land type scales in Northern Arkansas.

Fig. 7 .
Fig. 7. Simulated combined species density and basal area by age class at 1978, 2008, and 2128 in Northern Arkansas.

Fig. 8 .
Fig. 8. Simulated density, basal area, and biomass by species group, and carbon of total species at landscape and land type scales over 150 simulation years in Northern Arkansas.

Fig. 9 .
Fig. 9. Gingrich stocking charts (a) and Reineke density diagrams (b) showing mean stand trajectories from 1978 to 2128 (150 simulation years) at landscape and land type scales in Northern Arkansas.

Table 1 .
Seedlings establishment in stand initiation stage determined by species shade tolerance and growing space occupied (GSO) in the forest landscape model LANDIS PRO.

Table 2 .
Calibrated species life history parameters used in in the forest landscape model LANDIS PRO in Northern Arkansas.