The influence of land use and climate change on forest biomass and composition in Massachusetts, USA

. Land use and climate change have complex and interacting effects on naturally dynamic forest landscapes. To anticipate and adapt to these changes, it is necessary to understand their individual and aggregate impacts on forest growth and composition. We conducted a simulation experiment to evaluate regional forest change in Massachusetts, USA over the next 50 years (2010–2060). Our objective was to estimate, assuming a linear continuation of recent trends, the relative and interactive inﬂuence of continued growth and succession, climate change, forest conversion to developed uses, and timber harvest on live aboveground biomass (AGB) and tree species composition. We examined 20 years of land use records in relation to social and biophysical explanatory variables and used regression trees to create ‘‘probability-of-conversion’’ and ‘‘probability-of-harvest’’ zones. We incorporated this information into a spatially interactive forest landscape simulator to examine forest dynamics as they were affected by land use and climate change. We conducted simulations in a full-factorial design and found that continued forest growth and succession had the largest effect on AGB, increasing stores from 181.83 Tg to 309.56 Tg over 50 years. The increase varied from 49 % to 112 % depending on the ecoregion within the state. Compared to simulations with no climate or land use, forest conversion reduced gains in AGB by 23.18 Tg (or 18 % ) over 50 years. Timber harvests reduced gains in AGB by 5.23 Tg (4 % ). Climate change (temperature and precipitation) increased gains in AGB by 17.3 Tg (13.5 % ). Pinus strobus and Acer rubrum were ranked ﬁrst and second, respectively, in terms of total AGB throughout all simulations. Climate change reinforced the dominance of those two species. Timber harvest reduced Quercus rubra from 10.8 % to 9.4 % of total AGB, but otherwise had little effect on composition. Forest conversion was generally indiscriminate in terms of species removal. Under the naı¨ve assumption that future land use patterns will resemble the recent past, we conclude that continued forest growth and recovery will be the dominant mechanism driving forest dynamics over the next 50 years, and that while climate change may enhance growth rates, this will be more than offset by land use, primarily forest conversion to developed uses.


INTRODUCTION
Climate and land use changes are two major global ecological concerns with the potential to transform forest landscapes. The combined influence of these largescale anthropogenic forces is being superimposed onto naturally (and culturally) dynamic systems. Too often the ecological consequences of these forces are considered independently (Dale 1997). Scientists and policy makers require integrated analyses of multiple interacting processes in order to anticipate and adapt to future global change. We have conducted a series of landscape simulations that incorporate and project forward the current trends in forest growth and succession, land use, and climate change over the next 50 years  for the state of Massachusetts in the northeastern United States. Within the simulations, we monitored the relative influence of each of these factors on tree species composition and the stores of live aboveground biomass (AGB), which is the most dynamic and manageable pool of forest carbon (Fahey et al. 2010).
The current position of northeastern forests along their successional trajectory is largely attributable to the region's land use history ). The legacy of widespread agricultural clearing during European settlement (ca. 1600-1800) followed by old-field succession is thought to be a principle mechanism responsible for modern species structure and composition (Motzkin et al. 1996, Fuller et al. 1998, Cogbill et al. 2002, Hall et al. 2002) and rates of carbon storage and sequestration (Caspersen et al. 2000, Houghton 2003, Albani et al. 2006. The modern forest landscape is dominated by Manuscript  generalist and early-to mid-successional tree species such as red maple (Acer rubrum), white pine (Pinus strobus), and red oak (Quercus rubra), with lesser amounts of late-successional, longer-lived species such as beech (Fagus grandifolia), which dominated the precolonial forest. (Fuller et al. 1998, Cogbill et al. 2002, Hall et al. 2002. Although a diversity of conditions certainly exists, the region's typical forest is relatively young (80-120 years), aggrading, and even aged (Kelty et al. 2003). Continued forest succession, in the absence of disturbance, will drive changes in structure (larger trees) and composition as short-lived species are replaced by longer-lived and shade-tolerant species. Succession is also anticipated to be a dominant mechanism driving carbon sequestration over the coming decades (Albani et al. 2006, Keeton et al., in press). However, the scarcity of remaining old-growth forests to use as reference conditions makes predictions regarding future stand development uncertain. Conventional wisdom and ecological theory, largely born out of early forest simulation (gap) models, suggest that biomass will peak when forests are younger than 200 years and that ''steady state'' biomass dynamics will dominate thereafter (Bormann and Likens 1979). However, several recent studies of remnant old-growth stands suggest that maximum biomass will not be reached until much later, and that stands .350 years old have not reached equilibrium (Luyssaert et al. 2008, Lichstein et al. 2009). This suggests that forest growth and succession in the Northeast could continue to serve as a carbon sink for many decades into the future. Understanding the degree to which future forest dynamics in the Northeast will be attributable to continued growth and succession is complicated by anthropogenic climate change. Alterations in patterns of precipitation and growing season can interact with forest succession and affect species composition and biomass dynamics. Observed changes in climate over the past century in the Northeast have largely mirrored global trends, including a 0.88C increase in air surface temperature and more variable precipitation patterns (Hayhoe et al. 2007). Climate change projections for the region suggest that the temperature will continue to increase by 2.1-5.38C over the next century, depending on the level of increase in atmospheric greenhouse gas concentrations. Precipitation is projected to increase by as much as 30%, mostly in the winter months ). In the longer term (.100 years) climate change in the region will have a profound influence on tree migration, reproduction, and establishment success Prasad 1998, Mohan et al. 2009). Most tree species will shift northward or succumb to climate stress (including indirect effects such as changing disturbance regimes), increased competition, or other pressures. These studies notwithstanding, the full consequences of climate change and its indirect impacts on forests remain far from certain. The interacting effects of climate-related issues such as CO 2 fertilization, drought stress, and proliferation of forest pests, pathogens, and invasive species confound reliable predictions (Campbell et al. 2009, Dukes et al. 2009). In the shorter term (,100 years), most tree species are expected to persist in their current range even as their optimal climate envelopes shift northward. Beyond persistence, forests in the Northeast may experience increased growth rates and carbon sequestration owing to warmer temperatures, a longer growing season, and potential CO 2 -driven increases in photosynthesis and water-use efficiency (Ollinger et al. 2008). An optimistic perspective is that climate-related gains in carbon sequestration will act as a significant negative feedback in the climate system, with the potential to offset some of hazards posed by accelerating carbon emissions.
Any climate-related increase in carbon sequestration rates must be weighed against reductions in sequestration and storage capacity owing to changing land use patterns and intensity. Throughout much of the late 19th and all of the 20th century, the dominant land cover transition in the eastern United States was from agriculture to forest, which resulted in a substantial gain of forest area (Ramankutty and Foley 1999). While agricultural land cover continues to decline throughout the region, land cover transitions to developed uses now override reforestation and culminate in a net loss of forest cover (Drummond and Loveland 2010). Indeed, exurban forest conversion to low-density housing (6-25 homes/km 2 ) is the fastest growing form of land cover change in the United States. Between 1950 and 2000, the portion of rural low-density housing area increased from 5% to 25%, with some of the most rapid rates of forest conversion to developed uses found in Northeast forests (Brown et al. 2005). After more than a century of increasing forest area, each New England state is now losing forest cover (Foster et al. 2010). The effect of land use change on forest composition and carbon stores is regulated not only by its total footprint but also by its spatial distribution. For example, fragmentation of forest landscapes may reduce the amount of forest and alter forest patch configuration, which, in turn, reduces the density of mature, reproducing individuals that can regenerate (Iverson et al. 2004). In addition, forest losses in productive regions have a greater impact on longterm carbon dynamics than the equivalent area lost in less productive regions. Spatial patterns of land cover change, such as deforestation, are driven by socioeconomic as well as biophysical factors (Lambin et al. 2001). In the Northeast, like most industrialized regions, the rate and pattern of forest conversion to developed uses is often associated with economic growth and with the distance to urban centers and other amenities (Schneider and Pontius 2001, McDonald et al. 2006, Theobald 2010.
Less conspicuous land uses, such as timber harvest, also affect tree composition and forest carbon dynamics. Timber harvest intensities vary widely, but typically have a lesser impact on forest ecosystems than the ''hard deforestation'' associated with conversion to developed uses. Nonetheless, even partial harvest can affect species composition and ecosystem carbon budgets over time (Nunery and Keeton 2010), as well as the conditions for regeneration and successional change. A lack of highquality spatial data sets describing regional harvest activities often hinders efforts to assess the impact of timber harvest regimes on forest conditions. Although regional inventory data confirm that timber harvest is occurring across the landscape, specific characteristics and the spatial distribution of harvest activities are typically unavailable. Data are especially hard to obtain for private lands, which dominate forest ownership in the Northeast. Fortunately, a unique database exists for Massachusetts that describes the spatial pattern and volume of tree species removed for every cutting event conducted in the state over a 20-year period (1984McDonald et al. 2006). From this database we know that temporal trends in harvest have been consistent over this period, with total harvest volume typically ;450 000 m 3 /yr. Red oak and white pine are consistently the most harvested species. Spatially, the percentage of forests harvested varies widely across the state, ranging from 0.01% to 1.5% of forests annually, depending on the region. Harvest activity in Massachusetts, as in other regions dominated by private ownership, is most intense in rural areas and is generally negatively correlated with urbanization (Wear et al. 1999, Munn et al. 2002, McDonald et al. 2006. The effects of succession, climate change, urbanization, and timber harvest on forest ecosystems are complex and interacting. Scientists and policy makers need to understand the probable impacts of these forcings in combination. As an example, the 10-state cap-and-trade Regional Greenhouse Gas Initiative (RGGI), of which Massachusetts is a signatory, requires states to minimize their emissions by 10% by the year 2018 (RGGI available online). 6 To help achieve this goal, carbon sequestered in forests can be calculated and then sold to greenhouse gas polluters in order to meet 3.3% of their compliance obligation. For now, only afforestation programs are allowed; however, there may be opportunities to include reduced deforestation in the future if there are reliable ways to account for it. Analyses that integrate climate and land use change are an important step in that direction.
Spatially explicit landscape simulation models offer a tractable way to integrate multiple interacting and stochastic processes , and are particularly useful for examining scenarios of forest growth and land use across complex, multi-owner regions (Thompson et al. 2006). We used LANDIS-II ), a well-established forest landscape succession and disturbance model, to simulate regional forest dynamics as they are affected by climate change, timber harvest, and conversion to non-forest over the next 50 years. In this paper, our objective was to estimate, based on recent trends, the relative influence of continued growth and succession, climate change, forest conversion to developed uses (hereafter ''conversion''), and timber harvest on forest composition and aboveground forest biomass (and, therefore, carbon) in Massachusetts. In future analyses, we plan to use this parameterization of the model as a basis to explore alternative scenarios of landscape change, such as scenarios incorporating biomass energy, invasive insect defoliators, and forest conservation.

Study area
We examined forest dynamics throughout the state of Massachusetts, in the northeastern United States (69.9-73.58 E, 41.3-42.98 N; see Plate 1). Simulations were restricted to the 12 000 km 2 classified as forests within the land use database derived from the Resource Mapping Project at the University of Massachusetts, maintained by the Massachusetts Geographic Information System (available online). 7 Land-cover data were manually classified from 1:25 000 color aerial photography for the years 1985 and 1999. We focused on forest conversion to non-forest but acknowledge that some non-forest land may convert to forest in the future; however, this is, and will likely continue to be, a rare land cover transition and it was not considered in this analysis. Within our study, we utilized the U.S. Environmental Protection Agency's Level III ecoregions to delineate areas of distinct climatic and edaphic conditions (EPA 2010). There are 13 Level III ecoregions within the state (Fig. 1). The western ecoregions are mostly mountainous and forested, with granitic and metamorphic bedrock, with the exception of the Marbled Valley and Vermont Piedmont, portions of which are underlain by calcareous bedrock. Forest cover is typically lower (and urban cover typically higher) in the eastern ecoregions (especially Boston Basin and Cape Cod and Islands), where the topography is gentler. Within our simulations, we treated the forest environment within an ecoregion as homogeneous in terms of soil water-holding capacity and climatic conditions now in and the future (Table 1).

LANDIS-II model parameterization
We simulated forest and land use dynamics using the spatially interactive landscape model, LANDIS-II . LANDIS-II is designed for simulating forest dynamics over mesoscales (10 4 -10 7 ha), including establishment, competition, growth, decomposition, and biomass accumulation, while also integrating multiple disturbances such as wind, timber harvest, and, as shown for the first time here, forest conversion. Because LANDIS-II simulates forest dynamics over large scales, it does not track individual trees; instead, it utilizes species 3 age cohorts to achieve computational tractability ). Within LANDIS-II, forests are represented as a grid of interacting sites (cells). Each site is assigned to an ecoregion within which climate and soils are assumed to be homogenous. Within each site, the live aboveground biomass of each species 3 age cohort is tracked as it is influenced by growth, senescence, and various types and intensities of disturbance. Each site can serve as a seed source for mature species located on that site, with the probability of seed dispersal declining exponentially with distance from the site (Ward et al. 2005).
We modeled the 25 most abundant tree species as determined by stem counts in the U.S. Forest Service's Forest Inventory and Analysis (FIA) database . Species attributes such as shade tolerance, seeding distance, and sprouting ability were determined from the literature (see Table 2 for values and citations). The duration required for ANPP (aboveground annual net primary productivity) to reach maximum potential varied with species' shade tolerance; less shade-tolerant species achieved maximum ANPP faster, as indicated by the ''ANPP-shape'' parameter (Scheller and Mladenoff 2004). Within LANDIS-II, cohort biomass is a function of each species' maximum ANPP (maxANPP) and its maximum achievable live aboveground biomass (max-Biomass), the cohort's age, and inter-and intraspecific competition (Scheller and Mladenoff 2004). MaxANPP and maxBiomass are input parameters. We calculated maxANPP using the PnET-II generalized ecosystem process model (Aber et al. 1995) in a manner similar to that of Gustafson et al. (2009), Ravenscroft et al. (2010, and Scheller and Mladenoff (2004). The PnET family of models have been widely applied and validated in the northeastern United States (Aber et al. 1997, Ollinger et al. 1998, Ollinger and Smith 2005. To parameterize PnET-II, we obtained species specific estimates of foliar nitrogen concentrations and specific leaf mass for each species from the literature (see Table 2 for citations) and from the Northeastern Ecosystems Research Cooperative, Foliar Chemistry Database (available online). 8 All other inputs into PnET-II were derived from the generic conifer and northern hardwood values presented in Ollinger and Smith  . For the purposes of our simulations, climate and soils are assumed to be homogeneous within an ecoregion. Abbreviations are defined in Table 1.   2005). Estimates for the maximum achievable biomass for each species were determined from old-growth stands reported in the literature (see Table 2 for values and citations). Maximum biomass estimates were impossible to obtain for each ecoregion; therefore, relying on the known relationship between ANPP and maximum biomass (Keeling and Phillips 2007), we assigned the maximum reported biomass for each species to the ecoregion with the highest ANPP for that species and then scaled the other ecoregions downward proportionately (Table 2). Similarly, because maximum achievable biomass under future climates is unknown, we assumed that the relationship between ANPP and maximum biomass held and scaled linearly into the future. The initial forest condition was based on a spatial imputation of 591 field plots conducted by the FIA program using the gradient nearest neighbor method (B. T. Wilson and A. J. Lister, unpublished manuscript). The native resolution of the imputation map was 250 m, but we resampled to 100 m (using the nearest neighbor method) to better match the scale of the land use processes. We subset the imputation map to the area classified as forest in the MassGIS 1999 land use database. The FIA field plots represented within the imputation map only contain age information for one dominant tree ). To establish the species 3 age cohort information required to parameterize LANDIS-II, we estimated the age of all other trees based on their height and the site class using equations in Carmean et al. (1989). We adjusted the site class for each species based on conversion factors published by Dixon and Keyser (2008). Once ages were estimated for all trees, they were binned into 10-year species 3 age class cohorts for use in the model. At the onset of each simulation, LANDIS-II goes through a ''spin-up'' phase to calculate the starting biomass for each cell, wherein each species 3 age cohort is ''grown'' from establishment to its age at the beginning of the simulation (year 0). To evaluate our parameterization of LANDIS-II and its ability to simulate forest biomass dynamics in Massachusetts, we compared the biomass estimates as measured in the FIA field data to the LANDIS-II initial estimates (at year 0) using a Pearson's correlation coefficient and the root mean squared error (RMSE).
The probability that a species will establish on a given site (P est ) is a model input that varies spatially by ecoregion and temporally by changing climatic conditions. At a given site, simulated establishment is affected by light conditions, availability of seed or sprouting, and Notes: For shade tolerance, 1 is least tolerant, and 5 is most tolerant. Values for maxANPP and maxBiomass are given as means with range in parentheses. Maximum annual net primary productivity (maxANPP) was calculated for each of the 13 ecoregions using the PnET-II ecosystem model. Maximum biomass values were obtained from the literature and then scaled proportionately with maxANPP values. See Methods for details. Superscript numbers in parentheses refer to literature sources where values were obtained: 1, NERC (2010); 2, Ollinger and Smith (2005) for average northern hardwood and conifer values when values for foliar N were not available in NERC database or in the literature; 3. Smith and Martin (2001); 4, Reich et al. (1995); 5, Green et al. (2003); 6, Abrams and Kubiske (1990); 7, Richardson et al. (2001); 8, Woodwell (1974); 9, Kloeppel et al. (1993); 10, Lichstein et al. (2009); 11, maxBiomass was set to 250 Mg/ha when an estimate was unavailable in the literature; 12, Whitney (1993); 13, Pastor et al. (1984); 14, Vose and Swank (1993). the P est input parameter. Following methods described in detail within Scheller et al. (2005) and Gustafson et al. (2009), we used a simulation approach to develop P est . In brief, 300 years of weather conditions were stochastically simulated using the climate data, which were then used to develop the establishment modifiers (temperature, soil moisture, and minimum January temperature) developed by Pastor and Post (1985). These data, in combination with the species vital attributes, drought tolerance, cold tolerance, and soil nitrogen tolerance, were used to estimate P est for each species in each ecoregion and climatic condition.
We included small-scale wind disturbance within the simulations but, because it was not among the processes that we were explicitly considering in this study, we took a simple approach to its parameterization. Using the LANDIS wind disturbance module (Scheller and Mladenoff 2004), we set the wind rotation period to 350 years for inland ecoregions and 200 years for coastal ecoregions. The mortality probability associated with a wind event increased monotonically as a species cohort approached its maximum longevity. The average size of wind events was set to 0.1 ha, with few larger wind events allowed.
We performed a sensitivity analysis of nine key parameters within the model parameterization. We varied their values by 10%, or by one unit in the case of ordinal categorical variables, and assessed the corresponding impact on the total biomass at year 0 and year 50. This approach to sensitivity analysis is useful for gauging the relative importance of parameters and for assessing whether the state variables depend linearly on input parameters (Drechsler 1998).

Climate change
LANDIS-II was explicitly designed to simulate climate change effects on forested landscapes (Scheller and Mladenoff 2008). Within the LANDIS-II/PnET-II framework previously described, climate variables affect maxANPP, maxBiomass, and P est . We obtained all climate data (minimum temperature, maximum temperature, and precipitation) from the Northeast Climate Impacts Assessment Group (available online). 9 For scenarios without climate change, we used their average data from the period 1960-1999. For the scenarios that include climate change, we used projections for the period spanning 2010-2039 and 2040-2069 and divided the observed changes into decadal increments (Fig. 2). The data are from an average of three general circulation models, all portraying the Intergovernmental Panel on Climate Change (IPCC) B1 emissions scenario: U.S. National Oceanographic and Atmospheric Administration/Geophysical Fluid Dynamics Laboratory CM2.1 (Delworth et al. 2006 (Washington et al. 2000). The model average was statistically downscaled to one-eighth FIG. 2. Climate values used to estimate maximum annual net primary productivity (maxANPP), maximum biomass values (maxBiomass), and the probability that a species will establish on a given site (P est ). Dark lines represent the average values across the 13 Massachusetts ecoregions, and lighter shading depicts an envelope spanning from the minimum to maximum ecoregions. Climate data are from the Northeast Climate Impacts Assessment Group . See Methods for details. 9 hwww.northeastclimatedata.orgi degree by Hayhoe et al. (2008). A detailed comparison of predicted vs. observed values is given by Hayhoe et al. (2007). The IPCC B1 emissions scenario is among the most optimistic in terms of the degree of future warming (Nakicenovic 2000).

Land use
We made no attempt to predict the exact locations of future timber harvests or forest conversion. Instead, we created ''probability of harvest'' and ''probability of conversion'' zones, based on the past land use patterns in relation to several social and biophysical variables, which then dictated the spatial distribution of land use. To do this, we took a random sample of 10 000 points within the region classified by MassGIS as forest in 1985. In a GIS, we extracted information about each point that was described, whether it had since been converted to developed uses (based on the 1999 MassGIS Land Use Layer), or had been subject to timber harvest (based on state harvest permits during 1984(McDonald et al. 2006). We then identified a suite of potential predictor variables that have been associated with forest conversion and timber harvest in other studies conducted in the region (Table 3; see Schneider and Pontius 2001, Tyrrell et al. 2004, McDonald et al. 2006. In a GIS we extracted the values for all of the predictor variables at the location of the sample points. We used regression tree analysis (RTA), fitting one tree with forest conversion as the response variable and one tree with timber harvest as the response, to identify the relationships between the land use activities and the predictor variables. RTA is a nonparametric technique that recursively partitions a data set into subsets that are increasingly homogeneous with regard to the response (De'ath and Fabricius 2000). We used an implementation of RTA, called conditional inference trees, that is available within the PARTY library (Hothorn et al. 2009) within the R statistical language (R Development Core Team 2006). Conditional inference trees establish partitions based on the lowest statistically significant P value that is obtainable across all levels of all predictor variables, as determined from a Monte Carlo randomization test with 9999 permutations; this technique minimizes bias and prevents over-fitting and the need for pruning (Hothorn et al. 2006). We used the probability at terminal nodes of the regression trees, weighted by the proportion of area that they represented within the landscape, to define ''probability of harvest'' and ''probability of conversion'' zones throughout the study area.
We used the Massachusetts harvest permit records (McDonald et al. 2006) and an understanding of private landowner tendencies (Belin et al. 2005, Finley and to develop four different timber harvest prescriptions ( Table 4) that emulate the recent harvest regime in terms of their annual extent, size of harvest, and the species and age classes removed. State records indicate that the annual harvested area has ranged from 8500 to 14 000 ha/yr with no temporal trend. Accordingly, we simulated 10 500 ha of timber harvest per year, distributed across the state according to the probability zones we have described. Harvesting focused on older cohorts of economically valuable species (i.e., Pinus stobus, Quercus sp., Tsuga canadensis, and Acer saccharum).
LANDIS-II had not previously been used to simulate forest conversion. To accommodate this, we modified the existing timber harvest module (Gustafson et al. 2000) such that when a site is identified for conversion, a user-specified amount of forest biomass is removed across all species on the site, future establishment is prevented, and the maximum achievable biomass (maxBiomass) on that site is permanently reduced. For example, assume that a ''Small development'' prescription (0.25 ha) is scheduled to occur on a site at year 10, at which time the site contains 100 Mg of biomass and has a maximum achievable biomass of 200 Mg; at that point in the simulation, 25 Mg of biomass is removed, no new forests may establish on that site for the rest of the simulation, and the maximum achievable biomass is permanently reduced by 25% to 150 Mg. This approach allowed us to emulate patterns of rural forest conversion at scales smaller than our 1-ha grain size and it allowed continued, albeit reduced, biomass accumulation on sites that had experienced some level of conversion.
To estimate the annual amount of forest converted to developed uses, we divided the total number of hectares that were classified as ''forest '' in 1985 and ''developed'' in 1999 (total ¼ 68 432 ha) by the number of years between the interpretations (total ¼ 14 years), which resulted in 4888 ha/yr. This level of annual forest conversion was spatially distributed according to the RTA-derived probability zones we have described. Of course, forestland that is converted to developed uses often retains significant residual tree cover (Nowak and Crane 2002). To emulate the past level of residual tree cover within converted sites, we randomly selected 300 areas that had been converted to developed uses in the period between 1985 and 1999 and examined the proportion of residual tree cover within them using GoogleEarth. Within each polygon, we randomly placed 10 points and counted the proportion that overlayed tree canopies. We used the resulting distribution to help set the proportion of forest biomass removed within the land use footprint (Table 4).

The simulation experiment
To estimate the relative influence of forest succession, climate change, forest conversion to developed uses, and timber harvest on species composition and biomass, we conducted a series of simulations in a full-factorial design (where climate change, conversion, and harvests were incorporated as ''treatments''). This required eight different scenarios. A power analysis indicated that replicating each scenario five times was sufficient to capture the between-run variability. The low number of required replicates reflects the fact that the many stochastic components operating within LANDIS-II stabilize to their average when measured at the scale of the entire study area. We described and analyzed simulation results at years 25 and 50.
We used a three-way factorial ANOVA to assess the relative influence of treatments on total aboveground biomass, including the main effects and interactions of climate change, forest conversion, and timber harvest. We also summarized changes in biomass 3 species rank abundance and, to assess the relative influence of treatments on tree species composition, we used nonparametric multivariate analysis of variance (i.e., PerMANOVA) using the ''adonis'' function and a Bray Curtis distance matrix to describe community composition within the Vegan Community Ecology Package (Oksanen et al. 2009) in the R statistical language (R Development Core Team 2006). This technique is wellsuited for partitioning distance matrices among multiple sources of variation and fitting linear models to distance matrices, particularly when multivariate normality cannot be assured and there are more species than replicates (Legendre and Anderson 1999). The method uses a permutation test and pseudo-F ratios to estimate Notes: Development and harvest categories include (in parentheses) the percentage of the total area affected. Forest conversion prescriptions describe the size and intensity of forest removal simulated in LANDIS-II to emulate recent trends in forest conversion to developed uses. Timber harvest prescriptions describe the species and ages of forest cohorts removal and the size of the harvest units simulated in LANDIS-II to emulate the recent trend in state timber harvests. Each year the total area affected by all the forest conversion prescriptions is set to 4888 ha. Of that total, 50% (2444 ha/yr) is ''developed'' using the ''Small development'' prescription (simulating small house lots). Likewise, 25% of the total (or 1222 ha/yr) is affected by the ''Medium development'' prescription, and so forth. For the right-hand ''harvest'' regime, with 10 500 ha/yr affected, 33% harvest is equivalent to 3500 ha/yr. Species codes are defined in Table 2.
For any stand to be eligible for timber harvest, at least 40% of the stand must contain harvestable trees. P values. We performed separate PerMANOVA tests on the total AGB by species data and on the relative proportions of species' AGB.

Initial forest condition and model sensitivity
At the scale of the entire state, the LANDIS-II spinup of the initial forest condition resulted in an estimate of 181.83 Tg of AGB, which was within 15% of the total biomass estimated from the imputation of the FIA data conducted by B. T. Wilson and A. J. Lister (unpublished manuscript). At the plot scale, the Pearson's correlation between predicted (LANDIS-II) and observed (FIA field plots) AGB was 0.69 with a root mean square error of 44 Mg/ha (Fig. 3). There were no apparent biases associated with different dominant tree species within a site.
Overall, our parameterization of LANDIS-II was not overly sensitive to any of the individual parameters tested; that is, the percentage change in AGB associated with a 10% change (or one unit change for categorical variables) in a given input parameter was typically ,10% (Table 5). The ANPP shape parameter that controls the time that it takes a new cohort to achieve maxANPP was the most influential input parameter during the ''spin-up'' phase of the model (i.e., AGB at year 0). A 10% change in the ANPP shape parameter was associated with a 12-13% change in initial AGB. However, this effect was diminished by simulation year 50, when a 10% change was associated with a 4.5-6.5% change. Maximum ANPP was also a relatively influential parameter, with a 10% change resulting in a 7.6% change in AGB at year 0 and a 3.5-4.2% change in AGB at year 50. A 10% change in the maximum achievable biomass parameter had a relatively small effect at year 0 (2.2-2.5%), but by year 50 it was among the more influential parameters (5.8-6.5%). Adjustments of 10% (or one unit) in all other parameters resulted in ,2% changes in AGB at year 0 or year 50.

Land use
The RTA describing factors associated with forest conversion to developed uses identified six statistically significant (P , 0.001) partitions using six different predictor variables, which resulted in a tree with seven terminal nodes (Fig. 4) Table 2. this was the best predictor variable that we evaluated. The distance from forestland to non-forestland was the next most important predictor. The terminal node in the regression tree with the highest probability of conversion included those forests that were not protected, were close to existing built areas (,283 m), on shallow slopes (,7.9%), and in census blocks with comparatively high human population densities (.460 people/km 2 ). Forests meeting these criteria were predominantly in the eastern, coastal part of the state, with a lesser component in the urbanized Connecticut River Valley. Forests with the lowest probability of conversion (aside from those under permanent protection) were those characterized by greater distances to non-forest areas, with low housing densities (,36.9 houses/km 2 ) and low road densities (,1.3 km/km 2 ). The aerial photo interpretation of residual tree canopy cover within sites recently converted from forest to developed uses in the land cover map resulted in a canopy cover of 55% 6 22% (mean 6 SD; Fig. 5). There was no significant relationship between patch size and tree cover (P ¼ 0.214, R 2 0.04). We used this information to set the intensity of forest removal within the five forest conversion prescriptions described in Table 4.
The RTA describing factors associated with past timber harvests identified seven significant partitions (P FIG. 4. Regression tree used to create probability (color-coded percentages) of future forest conversion zones used within the LANDIS-II simulations. The regression tree model used the pattern of forest conversion observed in the period spanning 1985-1999 in relation to the predictor variables described in Table 3. Terminal nodes of regression trees were scaled by their area on the landscape and were used to spatially allocate forest conversion within simulations. Units for the predictor variables are given in Table 3.
, 0.001) using five predictor variables, which resulted in eight terminal nodes (Fig. 6). Road density was identified as the best predictor variable. Areas with a low density of roads (,2.17 km/km 2 ) generally had higher rates of timber harvest. Within this lower road density group, forestland within census blocks where the median home values were less than $171900 were associated with the highest probability of timber harvest. Most forestland meeting these criteria was west of the Connecticut River Valley ecoregion. Forestlands with higher road density (,2.17 km/km 2 ) and shallow slopes (,6.2%) or high median home values (.$187000) were associated with the lowest rates of harvest and were generally found in the eastern one-third of the state. Overall, the probability of timber harvest was strongly and negatively correlated with the probability of forest conversion.

Changes in AGB and species composition
Continued forest growth and succession resulted in net positive changes in total AGB throughout all simulations (Fig. 7). Excluding any of the three treatments (i.e., climate, timber harvest, and forest conversion), total AGB increased from 181.83 Tg to 259.56 Tg at year 25 (a 42% increase) and to 309.56 Tg at year 50 (a 70% increase) (Fig. 7). Of the three treatments, forest conversion to other uses had the largest effect on AGB (Table 6). Compared to the no-treatment run, forest conversion reduced AGB by 9.39 Tg at year 25 and by 23.18 Tg at year 50. Timber harvest also reduced total AGB when compared to the ''growth only'' simulation, albeit by a much lesser amount. By year 25, Timber harvests reduced total AGB by 3.89 Tg by year 25, and by 5.23 Tg by year 50. Climate change had a positive effect on growth, increasing total AGB by 5.33 Tg at year 25 and 17.3 Tg at year 50. The factorial ANOVA indicated that, at year 25, all of the main effects were significant (P , 0.001), but there were no significant interactions among treatments. By year 50, the main effects remained significant and there were also small, but significant, interactions between climate and the two land use variables. Differences in AGB that were attributable to the treatments strongly reflected the spatial distribution of land use. The largest differences between the ''growth only'' runs and ''current trends'' runs (which included climate change, conversion, and harvest) were within the Boston Basin, Bristol Lowlands, and Connecticut River Valley, respectively (Fig. 8).
Although the treatments had significant effects on total AGB, their effects on relative composition were comparatively minor. P. strobus and A. rubrum were ranked first and second, respectively, in terms of their total AGB throughout all simulations (Table 7). Climate change reinforced the dominance of those two species. The PerMANOVA of treatment effects on total biomass at year 25 and year 50 indicated that forest composition was most significantly impacted by forest conversion, followed by climate (Table 8). The PerMANOVA of treatment effects on relative species composition identified no significant differences (Table 8). Forest conversion was generally indiscriminate in terms of species removal. Timber harvesting reduced Q. rubra from 10.8% to 9.4% of total AGB, but otherwise had little effect on composition (Table 7).

Biomass and species composition
In these simulations, the largest changes in AGB (and therefore carbon) arose from continued stand growth and succession, irrespective of any climate or land use changes. This finding reflects the fact that the average forest carbon density present in the landscape is well below what is seen in old-growth forests. This process of continued AGB accrual is operationalized within the model vis-a`-vis the initial forest biomass condition, which is based on forest inventory plots, and the maximum biomass parameters, which are based on empirical estimates taken from a sparse sample of oldgrowth forests, then adjusted downward based on estimated productivity of the ecoregion. Our analysis suggests that sustained forest recovery, owing to the legacy of agricultural abandonment, reductions in widespread forest harvesting, and regrowth following the 1938 hurricane and associated timber salvage, will continue to be the most important mechanism affecting forest carbon dynamics. Indeed, even assuming a continuation of the current trends in land use and FIG. 5. Distribution of tree cover within 300 randomly selected areas identified by MassGIS (a Massachusetts state government agency that disseminates geographic data; see hhttp://mass.gov/mgisi) as having been converted from ''forest'' to ''developed uses'' between 1985 and 1999. Tree cover was measured by counting the proportion of 10 randomly positioned points overlaid on a tree canopy within aerial imagery in GoogleEarth.
climate change, the simulations predicted a 65% increase in AGB over the coming 50 years. This finding, if born out, has significant implications for climate change mitigation and climate policy.
Of course, our modeling framework is simplistic as compared to actual forest landscape dynamics, and the capacity of the northeastern forest carbon sink remains a subject of considerable uncertainty and interest. One important counter perspective comes from long-term studies at the Hubbard Brook Experimental Forest in central New Hampshire, where researchers have found that forests recovering from intensive harvest, dating to FIG. 6. Regression tree used to create probability (color-coded percentages) of future timber harvest zones used within the LANDIS-II simulations. The regression tree model used the pattern of forest harvest observed in the period spanning 1985-2004 in relation to the predictor variables described in Table 3. Terminal nodes of regression trees were scaled by their area on the landscape and were used to spatially allocate timber harvest within the simulations. Units for the predictor variables are given in Table 3. the early 20th century, have stopped accumulating biomass at ;70-80 years, having accrued just 95-100 Mg/ha (Fahey et al. 2005), which is well before and well below the saturation point expected within our simulations. Notably, the northern hardwood forests at Hubbard Brook contain few biomass-dense species (such as Quercus spp.) and have biogeochemical limitations resulting from shallow soils and decades of acid rain (Campbell et al. 2007). A contrasting example comes from slightly older forests just 100 km to the south at the Harvard Forest in central Massachusetts, where the forests continue to accumulate biomass at a rate of 4-6 MgÁha À1 Áyr À1 , a rate that appears to be an accelerating (Foster et al. 2010; S. C. Wofsy and J. W. Munger, personal communication). The Harvard Forest plot is on a rather low-productivity site compared to the average productivity within the state of Massachusetts, which suggests that there is a potential through most of the study area for continued carbon storage over the coming decades. Other modeling and empirical studies support this perspective. For example, the length of the recovery period within our simulations is consistent with predic-  tions for the entire eastern United States made by the Ecosystem Demography (ED) model (Albani et al. 2006), which estimates maximum biomass as an emergent property of ecosystem processes, rather than a user input based on old-growth data. Expectations based on the ED model have biomass increasing at least until the end of the 21st century (Albani et al. 2006). There are also observational studies that have examined old forest and suggested a protracted biomass recovery period in the Northeast. For example, a study of remnant old-growth northern hardwood-conifer stands in New York showed that old-growth stands (250-400 years) contain .40% more aboveground biomass than mature stands (100-150 years) (Keeton et al., in press). Similarly, Brown et al. (1997) and Lichstein et al. (2009) have compared the upper end of the biomass distribution within the FIA database to values from old-growth stands and have concluded that there is potential for the landscape to accrue significant quantities of additional biomass. Because we did not simulate hurricanes, ice storms, or insect outbreaks, and only included a moderate amount of wind disturbance, actual biomass accrual over the next 50 years will probably lag behind our simulations. Nonetheless, given the range of supporting evidence, the pattern of forest biomass accrual owing to continued forest growth and succession that was estimated by these simulations seems reasonable.
After growth and succession, forest conversion to developed uses had the largest impact on AGB stores. By year 50, forest conversion reduced the total AGB by 23.18 Tg compared to the growth-only run. This effect was more than four times that of timber harvest, despite the fact that timber harvest occurred over twice the land area annually. This finding emphasizes the full costs of permanent forest conversion. Whereas timber harvest removes a portion of the forest biomass stored on a site, it does not necessarily affect the longer-term capacity for growth. In contrast, permanent forest conversion removes the stored biomass and the capacity for future growth. This is highlighted by the diminishing incremental effect from year 25 to year 50 for timber harvest (À3.9 to À5.2 Tg) vs. the increasing incremental effect for forest conversion (À9.4 to À23.2 Tg; Table 6). It is also worth noting that forest conversion occurred predominately in the eastern ecoregions, where forest productivity was the lowest, whereas timber harvest tended to be concentrated in the more productive ecoregions in the western part of the state. Had our simulations been aspatial or had we distributed land use randomly across the state, the effect of forest conversion on AGB stores would have been higher and the effect of harvest would have been lower.
The LANDIS-II simulations are consistent with others that suggest that future anthropogenic climate change, at least as it affects temperature and precipitation, will have a net positive impact on growth and biomass stores (Campbell et al. 2009). Total AGB at year 50 was 5.5% higher when climate change was simulated than when climate was held static. Despite this finding, the full effects of climate change remain far FIG. 8. Percentage of year-0 aboveground biomass at simulation year 50 by ecoregion for the ''Current Trends'' scenario, which included forest growth, climate, conversion, and timber harvest (top number) and for the ''Growth Only'' scenario, which included only forest growth (bottom number) of each pair. from certain. On the one hand, CO 2 fertilization, which we did not consider, could accelerate growth even above what we suggest here (Ollinger et al. 2008); on the other hand, increasing variability in environmental parameters such as sunlight, precipitation, and temperature could lead to lower rates of productivity that override the broadscale trends (Medvigy et al. 2010). Also, potential indirect effects of climate change, such as increased hurricane intensity, ice storm frequency and intensity, and the proliferation of insect pests, could reduce or reverse any climate-related increases in biomass storage. The largest sources of uncertainty, with regard to climate change effects on forest growth, are the emissions and climate predictions. For simplicity's sake, we only modeled one emission scenario: the IPCC B1. In this scenario, it is assumed that CO 2 concentrations are stabilized near 550 ppm (Nakicenovic 2000), which is a rather optimistic perspective, given the societal response to climate forecasts thus far. Although uncertainties seem to dominate our consideration of climate impacts on the forest biomass of Massachusetts, it seems likely, even given our somewhat simple portrayal, that any positive feedbacks associated with climate changerelated increases in carbon sequestration will not be sufficient to override the losses of forest biomass due to land use change.
With regard to tree species composition, our simulations suggested only minor changes in relative abundance, irrespective of the effects of climate or land use change. This finding largely reflects the difference between the longevity of the trees (typically . 200 years) relative to the length of the simulation (50 years). Many studies have demonstrated the important impacts that climate change will have on tree species distribu-  Notes: Values within parentheses are percentage of total AGB contributed by each species. Species codes are defined in Table 2. ''CC'' is climate change, ''SC'' is static climate, ''FC'' is conversion of forest to developed uses; ''Harv.'' is timber harvest. tions in the Northeast, but these typically have model horizons out 100 years or more (e.g., Iverson et al. 2004). In general, our simulations suggest that the landscape is very slowly shifting dominance to latesuccessional species (i.e., Tsuga canadensis, Fagus grandifolia). Also significant in our simulations was the impact of the timber harvest regime on Quercus rubra. At year 50, in the full current trends runs, Q. rubra makes up 8.6% of the forest AGB, but with harvesting removed Q. rubra increases to 10.1% (a difference of 4.9 Tg of AGB). Q. rubra is a valuable timber species and is consistently among the most sought after in Massachusetts (McDonald et al. 2006). The effect of timber harvest on its landscape abundance, while marginal, should be considered in light of the region-wide trend of diminishing recruitment of Quercus, as evidenced in the FIA database and elsewhere (e.g., Abrams 2003, McEwan et al., in press). Further study is needed to understand how (or if ) the timber harvest regime is exacerbating the other causes for regional oak decline. The full utility of having specieslevel interactions in this parameterization of LANDIS-II will come in our next stage of research, when we will incorporate scenarios of invasive pests and pathogens (i.e., the hemlock woolly adelgid, Adelges tsugae) that affect particular species.

Limitation and challenges of simulating coupled natural and human systems
Our simulations are not to be interpreted as predictions or forecasts of any kind. Coupled natural and human systems are characterized by complex interactions, feedbacks, and time lags, which tend to manifest as surprises, rendering true prediction impossible (Liu et al. 2007). Rather, the simulations represent the decomposition of one rather basic scenario that portrays a linear continuation of current trends in land use and climate change. ''Current trends'' or ''business as usual'' scenarios are popular and are useful as a type of null model or straw man from which we can build suites of alternative scenarios, which may then help to bound the range of plausible futures. The ''current trends'' represented in our simulations capture a period in history (1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999) when the economy of Massachusetts was growing quickly and commercial and residential development rates were comparatively high (DeNormandie 2009). As such, the rate of forest conversion within these simulations should be interpreted in context of the robust economic environment. If, for example, the ''current trends'' emulated the period between 2008 and 2010, when a national recession slowed rates of new building construction, there would be far less forest conversion. The rate and intensity of timber harvest, in contrast, has been relatively consistent over time and it appears to be less sensitive to the larger economic context (McDonald et al. 2006). Nevertheless, past stability could belie future trends in timber harvest if, say, the current initiatives to build biomass energy plants were to gain traction.
Our approach for modeling the spatial distribution of land use differs from most land use change simulations, which tend to focus on precise spatial allocations of land change, where the probability of change is dynamic throughout the simulations (e.g., Silva andClarke 2002, Verburg et al. 2002). These approaches explicitly acknowledge that patterns of land use are reactive to land use in the time steps before. We took a coarser approach and defined probability of land use zones that were static over the duration of the simulation and within which land use was allocated randomly. This approach, which does not permit creeping patterns of sprawl, resulted in densities of forest conversion within the high-probability zones that were probably unrealistically high by the end of the simulations. However, given that our objectives were outside any efforts in land use planning, this approach was an effective way to capture the broad spatial structure of land use. Moreover, efforts to model the precise spatial distribution of future land use are more often wrong than right (Pontius et al. 2008).
In addition to the uncertainties associated with the land use and climate change scenarios and the simplified approach to allocating land use, there are other limitations of our approach that must be acknowledged. Importantly, the parameterization of the LANDIS-II and PnET-II models are, by definition, simplifications of complex processes and they are limited by the data available in the literature. At their best, models such as these should be interpreted as formalizations of the current state of knowledge regarding several interacting processes. To reiterate, our results should not be interpreted as predictions or forecasts of the future; instead they highlight the relative importance of different processes and offer a platform for testing assumptions.

CONCLUSIONS
The continued growth of forests within Massachusetts, like much of the eastern United States, has a strong element of inertia that has been building since the era of agricultural abandonment (circa 1850-1900). Our simulations suggest that this century-old legacy of land use change will remain the dominant mechanism controlling forest biomass and tree compositional dynamics for at least the next 50 years. To be sure, the modern land use regime does affect the forest landscape, and these impacts reduce the amount of forest biomass and have some impact on tree species composition.
What is more, in the case of forest conversion, the loss is significant in terms of both the biomass removed and the loss of capacity to recover and grow in the future. However, the modern land use regime affects less than 1.5% of the forested landscape per year, and within that footprint much of the forest cover remains. As a result, modern land use pales in comparison to the inertia of forest growth. Our simulations suggest a similar story with regard to the influence of anthropogenic climate change on forest biomass and tree composition. It seems likely that climate change will have a positive effect on growth and carbon sequestrations rates if, as the climate models suggest, growing seasons lengthen and precipitation rates increase. However, the effect of climate change will be small relative to the background rate of growth that is attributable to the age and vigor of the forest. Although climate change is undoubtedly shifting the optimal establishment windows for tree species, for now and for the next 50 years, the existing forest will largely just continue to grow. Of course, over longer time frames or in the face of large disturbances (e.g., hurricanes or infestations), compositional shifts will be more apparent.