Climate Change Effects on Northern Great Lake (USA) Forests: A Case for Preserving Diversity

Under business as usual (BAU) management, stresses posed by climate change may exceed the ability of Great Lake forests to adapt. Temperature and precipitation projections in the Great Lakes region are expected to change forest tree species composition and productivity. It is unknown how a change in productivity and/or tree species diversity due to climate change will affect the relationship between diversity and productivity. We assessed how forests in two landscapes (i.e., northern lower Michigan and northeastern Minnesota, USA) would respond to climate change and explored the diversityproductivity relationship under climate change. In addition, we explored how tree species diversity varied across landscapes by soil type, climate, and management. We used a spatially dynamic forest ecosystem model, LANDIS-II, to simulate BAU forest management under three climate scenarios (current climate, low emissions, and high emissions) in each landscape. We found a positive relationship between diversity and productivity only under a high emissions future as productivity declined. Within landscapes, climate change simulations resulted in the highest diversity in the coolest climate regions and lowest diversity in the warmest climate region in Minnesota and Michigan, respectively. Simulated productivity declined in both landscapes under the high emissions climate scenario as species such as balsam fir (Abies balsamea) declined in abundance. In the Great Lakes region, a high emissions future may decrease forest productivity creating a more positive relationship between diversity and productivity. Maintaining a diversity of tree species may become increasingly important to maintain the adaptive capacity of forests.


INTRODUCTION
Over the next century, the pace of climate change is projected to exceed the ability of forests to naturally adapt (Soja et al. 2007, Iverson et al. 2008, Loarie et al. 2009, Kuparinen et al. 2010, Svenning and Sandel 2013).For example, continental declines in quaking aspen (Populus tremuloides) throughout North America are linked to climate change (Worrall et al. 2013).Generally, substantial shifts in the ranges of tree species are expected under climate change (Zhu et al. 2012).
However, a lag in species response may result in tree species extirpation, loss of forest structure and productivity, a decline in species diversity, and ultimately a decline of ecosystem resilience (Soja et al. 2007, Scheller et al. 2008, Frelich and Reich 2010, Fisichelli et al. 2012).The changes may exceed land-owners' ability to successfully manage forests with traditional forestry strategies.Within the northern Great Lakes region, vulnerability to climate change will likely depend on both latitude and proximity to the boreal-temperate transition zone (Fisichelli et al. 2013).More specifically, adaptation to climate change may depend on the potential mix of species capable of providing ecosystem services under changing environmental conditions.''Insurance'' species are those that may currently occur in low abundance however have a response capability (i.e., the ability to provide ecosystem services) beyond the range of environmental tolerances provided by the more dominant species (Naeem andLi 1997, Walker et al. 1999).
The northern Great Lakes region has and is expected to experience fast rates of climate change compared to continental edges where climate change is expected to be more buffered by oceans (IPCC 2007).Additionally, the upper Midwest represents an ecotone in the borealtemperate forest transition zone where many species are at the southern edge of their ranges (Curtis 1959).By late in the 21st century, mean winter temperatures in the Great Lakes region are projected to increase by 3-78C, and mean summer temperatures are projected to increase by 3-118C (Kling et al. 2003).In addition, longer growing seasons and greater variability in precipitation are also expected (IPCC 2007).Current greenhouse gas emission observations have been measured at or above the highest projected emissions scenario (A1FI) provided by the Intergovernmental Panel on Climate Change (IPCC 2007, Peters et al. 2012).These results give interest and urgency to understanding future forest dynamics in the Great Lakes region.The Global Circulation Model (GCM) developed by NOAA's Geophysical Fluid Dynamics Laboratory (GFDL) using the A1FI ('high') emissions scenario and future downscaled climate projections throughout the Great Lakes region indicate that temperature will increase in all seasons throughout the year (IPCC 2007, Stoner et al. 2012).In addition to projected future temperature increases, observed increases in temperature over the last century corroborate projections (Andresen et al. 2012).
Having diverse species with a broad range of environmental tolerance may build forest resistance (the ability to persist through disturbance with little change) and resilience (the ability to rebound following disturbance) to the effects of changing environmental conditions (Folke et al. 1996, Chapin et al. 1998).Some studies suggest that the relationship between diversity and productivity is more positively correlated in less productive sites compared to highly productive sites (Loreau et al. 2001, Whittaker et al. 2001, Paquette and Messier 2011).In productive temperate forests with deep, rich soils, competitive exclusion allows single species to dominate and increase overall productivity (Paquette and Messier 2011).In less productive boreal forests, positive diversity-productivity relationships have been found indicating less competitive exclusion (Lei et al. 2009, Paquette andMessier 2011).Nevertheless, even in highly productive forests, homogeneous forests may be vulnerable to the effects of a single host pathogen (e.g., bark beetle damage [Bentz et al. 2010]).If a loss of biodiversity is experienced, large impacts to productivity may result (Tilman et al. 2012).
The resilience of northern Great Lakes forests is, in part, generated by their natural disturbance regimes (Frelich 2002).The combination of harvesting, wind, and fire play a distinct role in facilitating the function of the extant major forest types.Within the northern Great Lakes forests, the most frequent and severe disturbances come from timber harvest by diverse land owners (Fig. 1) (Karamanski 1989, Heinselman 1996 ).Through a long history of harvesting in the region, these forests are now characterized by denser stands than their historic counterparts and are more homogeneous (Hanberry et al. 2012).In recent decades, an overall decline in harvest intensity and rate has been observed (Woodall et al. 2011), along with an increase in partial harvest systems (D'Amato et al. 2009).
Generally, wind and fire events in the Great Lakes region are low in frequency and small in size (Lorimer and White 2003, Schulte and Mladenoff 2005, White and Host 2008).Frelich's estimated annual percent area damaged by wind was less than 0.05% in Northeastern Minnesota, and in northern lower Michigan between 0.05% and 0.1%.White and Host (2008) estimated that blowdown events historically occurred on between 0.02% and 0.29% of the landscape annually in northeastern Minnesota.White and Host (2008) estimate percent area burned annually in Northeastern Minnesota to be between 0.1% and 0.6% of the landscape.In northern lower Michigan, Cleland et al. (2004) estimate the percent area burned annually to be between 0.09% and 0.13%.
Associations of tree species ('forest types' or 'species assemblages'), within the Great Lakes region include a mix of boreal and temperate forest types with specific threats, adaptation to disturbances, and vulnerability to climate change (Handler et al. in press a, Handler et al. in press b).Given the value of ecosystem services specific to species assemblages (e.g., habitat for the federally listed Kirtland warbler [Dendroica kirt-landii] in jack pine [Pinus banksiana] barrens [Walkinshaw 1983]), or differences in carbon storage and productivity between different forest types (Gheorghe et al. 2010), the persistence of unique forest types provides insurance for the production of future ecosystem services.Under changing conditions, a diverse suite of forest types may help to maintain or increase productivity (Hooper and Vitousek 1997).
Our objectives were to assess how climate change projections under business as usual (BAU) forest management may affect northern forests.Specifically, we addressed the following questions: (1) With the risk of declining productivity, how might climate change affect the relationship between species diversity and forest productivity?(2) Given the value of biodiversity, how might regional differences in climate change, soil productivity, and forest management affect tree species diversity in the Great Lakes region?(3) Because many ecosystems servicesincluding productivity-are dependent upon v www.esajournals.orgtree species and assemblages, we also asked how might specific tree species and forest types respond to climate change?

Study areas
Northeastern Minnesota and the northern lower peninsula of Michigan provide two examples of northern Great Lakes forested landscapes that represent a range of forest types and conditions.The Minnesota landscape is largely un-fragmented within the boreal-temperate ecotone (Curtis 1959).With the exception of the Boundary Waters Canoe Area (BWCA) and several protected state parks along Lake Superior, northeastern Minnesota is largely actively managed for timber (D'Amato et al. 2009).In contrast, northern lower Michigan is situated within the northern temperate forest.Northern lower Michigan forests include less actively managed forests but also fewer large tracts of protected contiguous forest land (Gustafson andLoehle 2006, The Conservation Biology Institute 2010).Both landscapes include a large proportion of private non-industrial land whose goals and management regimes vary (Potter-Witter 2005).We chose the landscape extents based on ecological boundaries comprised of 1.6 and 2.6 million hectares of forest land in Minnesota (White and Host 2000) and Michigan (Albert 1995), respectively (Fig. 1).

Experimental design
Our experimental design includes BAU management in our two landscapes under three climate scenarios (current climate, low emissions future, and high emissions future; see below).We used as consistent an approach as possible to parameterize the model across both the Minnesota and Michigan landscapes.For all simulations, we used a two-hectare resolution and a 150-year simulation horizon starting at year 2000.
We modeled forest change using the LANDIS-II (version 6.0) forest landscape model (Scheller et al. 2007).LANDIS-II is a widely used forest simulation model that simulates successional processes, including disturbance, seed dispersal, growth, and mortality.These processes are spatially dynamic and interact across the landscape.The landscape is made up of interconnect-ed cells within zones or 'ecoregions' of homogeneous climate and soil.Within each cell, regeneration, growth, and mortality occurs and trees are represented as species age-cohort combinations.The LANDIS-II framework is comprised of user-selected extensions that vary in process and complexity.We utilized the Biomass Succession extension (version 3.1) (Scheller and Mladenoff 2004), Biomass Harvest (version 2.1) (Gustafson et al. 2000), Base Fire (version 3.0) (He and Mladenoff 1999), and Base Wind (version 2.0) (Scheller and Mladenoff 2004).We ran these extensions at a 5-year time step.
Seed dispersal is an important process within LANDIS-II.A cell can act as a seed source if species-cohorts have grown to reproductive maturity.Seed dispersal probability declines exponentially from the edge of the source cell allowing 5% probability assigned to long distance seed dispersal and 95% probability assigned to short distance dispersal (Scheller et al. 2007).Short and long dispersal distances are species dependent (Burns and Honkala 1990).
We used the PnET-II forest ecosystem process model to calculate LANDIS-II species-specific input parameters (Xu et al. 2009), including maximum aboveground net primary productivity per year (ANPP max ), and probability of establishment (P est ) at every time step.These parameters are the primary mechanisms whereby climate change was incorporated into the simulations.We calculated P est using monthly maximum and minimum temperature, precipitation, and photosynthetic active radiation (PAR) within unique soil regions in each landscape.We calculated the Maximum Aboveground Biomass (AGB max ) based on the species-specific relationships observed between AGB max and ANPP max (Thompson et al. 2011).Other inputs to calculate PnET-II parameters included percent foliar nitrogen content, maximum foliar mass area, and soil water holding capacity (Table 1).
ANPP max , AGB max , and P est are input parameters to the Biomass Succession extension which grows, matures, and causes age-related mortality of species cohorts (Scheller and Mladenoff 2004).ANPP max determines the growth potential of above ground biomass moderated by age and competition (Scheller and Mladenoff 2004).P est determines the establishment of a new cohort, given a seed source (or seedling source in the case of planting) and adequate light, based on a climate envelope approach (Xu et al. 2009).At the beginning of each simulation, Biomass Succession ''spins-up'' initial biomass based on the simulated past growth of each species-age cohort, dependent on their age at the beginning of the simulation.Growth is limited by a species specific growth curve parameter which determines how fast actual ANPP achieves ANPP max .Maximum biomass of a species-cohort is further limited by competition within a cell as represented by potential growing space, minus space occupied by other species-cohorts.Age-related mortality is represented by an increasing decline in AGB close to the age of maximum lifespan by species.Following age-related mortality, a cohort's biomass is added to a detritus pool where a species-specific decay rate decomposes the biomass and where each species detritus becomes available to be altered by disturbance.

Ecoregion delineation
Within the landscape extents, we identified relatively homogeneous soil regions within climate regions to create unique ecoregions.In Minnesota, we used seven previously defined climate region boundaries (White andHost 2000, Ravenscroft et al. 2010).In Michigan, we used a cluster analysis to delineate nine relatively homogeneous climate regions from the PRISM observed climate database (Daly and Gibson 2002).The inputs for the cluster analysis included 30-year average monthly minimum temperatures for April, May, June, September, October, November, and average annual precipitation.We chose the fall and spring months to delineate temperature as they show large climate variability around thresholds for photosynthesis (i.e., Table 1.Tree species life history characteristics used in the LANDIS-II simulations.
Species simulated in the Michigan landscape only.à Estimated from Pinus banksiana.
§ Estimated from Populus tremuloides.} Species simulated in the Minnesota landscape only.
v www.esajournals.orgphotoperiod, freeze, and thaw temperatures).To delineate soil regions, we calculated soil water holding capacity (SWHC) values (cm of water within one meter of soil depth) in northern lower Michigan (SSURGO Soil Survey Staff 2011) and northeastern Minnesota (STATSGO 1994).The SSURGO database provided a preferable finer resolution however was unavailable within most of the northeastern Minnesota landscape.Within each landscape, we binned SWHC into 10 soil regions in ArcGIS (ESRI 2011).Finally, within each landscape, we nested our soil regions within climate regions to create unique ''ecoregions'' (Fig. 2).

Initial communities
Inputs to LANDIS-II include an initial conditions map with tree species-age cohorts assigned to every forested cell.For the Minnesota landscape, we utilized previously developed initial community data (Ravenscroft et al. 2010).For the Michigan initial community data, we imputed FIA plot data spatially to a combination of two forest cover maps: We leveraged the Forest Biomass Information System (FBIS) (Falkowski, unpublished manuscript) imputation map of stand age and Forest Inventory and Analysis (FIA) forest type (Falkowski et al. 2010).Likewise, we utilized an additional Michigan Department of Natural Resources forest classified map from the Integrated Forest Monitoring, Assessment, and Prescription Stage 1 map (IFMAP) (MDNRS 2001).Both maps were provided at 30-m resolution.To resample to two-hectare resolution and preserve the proportion of forested land within the landscape, we used a combination of the ''majority'' rule for contiguous regions (forested cells) and ''nearest neighbor'' rule for fragmented regions (non-forested cells) to resample the landscape within ArcGIS (ESRI 2011).As the IFMAP forest types were less aggregated than the FIA/FBIS forest types, we identified the one-to-many (IFMAP forest type to FBIS forest types) relationship with a look-up table.Assuming that the FIA forested plots from the northern lower peninsula of Michigan represent the distribution of forest types present, we adjusted the look-up table to best match percent of cells assigned to a given forest type with percent of the FIA plots with that forest type.Based on the final look-up table, we assigned an IFMAP forest type to every FIA forest type.Although the FBIS map provided more precise forest type detail, the accuracy of the data was less robust than the IFMAP forest type.For every cell with a matching FBIS and IFMAP forest type, we randomly imputed a FIA plot with the matching FBIS Forest type and FBIS stand age.For cells whose forest type values did not match between FBIS and IFMAP, we randomly imputed a FIA plot with the matching combination of: IFMAP forest type, FBIS stand age, and county from where the plot came.
To populate the initial communities map, we selected all the forested FIA plots from within the northern Lower Peninsula of Michigan and converted individual trees to species age cohorts.We used site index curves to estimate the age of each tree in our FIA plots (Eq.4.7.2.2 in Dixon and Keyser [2008] based on Carmean et al. [1989]).When multiple site index curves were provided for a species, we choose the curve coefficients with the closest geographic match to the northern Great Lake region.This method resulted in a better fit of FIA estimated biomass to LANDIS-Year 0 biomass compared to multiple regression equations using height, diameter, and crown class to predict tree age based on the assumption that stand age is equal to the average age of trees in the plot.Finally, we aggregated calculated tree ages to five-year age cohorts.

Climate data
We simulated three climate futures.We utilized greenhouse gas emission scenarios as representative of plausible futures.Emission scenarios are based on specific future changes in population, technology, and environmental priorities.We simulated current climate based on randomly selecting simulation years from 30 years of past observed climate   (Daly and Gibson 2002).We simulated a low emissions climate future based on the Parallel Climate Model (PCM) GCM (Washington et al. 2000) and the B1 emission scenario (IPCC 2007).The B1 emission scenario represents the lowest future increase in fossil fuel use described by IPCC.We simulated a high emission climate future based on the GFDL GCM (Delworth et al. 2006) and the A1FI emission scenario (IPCC 2007).The A1FI emission scenario represents the most fossil fuel intensive future described by IPCC and represents the closest emission scenario to current observed global emissions (Raupach et al. 2007).We chose these emission projections to bracket a range of uncertain climate futures.PCM climate projections are considered to have low sensitivity to emissions while GFDL projections are considered to be more moderately sensitive to emissions.We chose the combination of the B1 emission scenario coupled with the PCM global climate model and the A1FI emission scenario coupled with the GFDL global climate model to bracket a large range of plausible futures.These emission and GCM combination scenarios are also being used in other coordinated research efforts covering the region and will ensure consistency among model projections (Peters et al. 2013) The IPCC climate projections include climate futures to year 2100 (IPCC 2007).To simulate forest change to 2150, we interpolated climate variables (maximum and minimum temperature, precipitation, and PAR) for the remaining 50 years based on trends and variability in each climate variable within the original 100 years of data using the Amelia library in R (Honaker et al. 2011).Climate projections under both emission scenarios indicate increased temperature in all seasons over the next 100 years (Fig. 3).As expected, the lower latitude Michigan landscape and the higher emission climate scenario indicate higher temperatures than the Minnesota landscape and the low emission scenario, respectively.Precipitation is more variable; nevertheless a decrease in precipitation is projected under the high emission climate scenario in summer months after 2050 (Stoner et al. 2012).
In order to aggregate climate regions for analysis, we calculated average July temperatures for each climate region in each landscape from 30 years of historic observations   (Daly and Gibson 2002).We then grouped each of the 9 climate regions in the Michigan landscape and the 7 climate regions in the Minnesota landscape into 5 aggregated and ordered climate regions in each landscape based on the mean July temperature.For analysis of climate change, the aggregation and ordering assumes that the regions would stay ordered relative to each other (e.g., the warmest regions within a landscape will be the future warmest regions within a landscape).The warmest aggregated climate region in the Minnesota landscape was cooler than the coolest aggregated climate region in the Michigan landscape (Fig. 4).

Disturbance regimes
To estimate the wind and fire patch size and disturbance rotation period, we used historical assessments of disturbance in Minnesota (White and Host 2008) and Michigan (Whitney 1986, Zhang et al. 1999, Cleland et al. 2004).We assumed constant fire and wind disturbance regimes for each climate scenario.Compared to timber harvesting in the region, wind and fire disturbances represent minor disturbance regimes.We adjusted the mortality parameters of American beech (Fagus grandifolia) and ash (Fraxinus spp.) due to the current and expected pathogens beech bark disease and emerald ash borer, respectively.As these vectors primarily affect the older age cohorts for ash (Kashian and Witter 2011) and beech (Busby and Canham 2011), we adjusted the mortality shape parameter to allow the species to live to their maximum longevity but reduced their maximum biomass at an earlier age.We calibrated these values with the LANDIS-II-Site v2.3 tool (B.Miranda, personal communication) which is designed to simulate succession on a single site using parameters from We simulated harvesting and management with the Biomass Harvest extension (Gustafson et al. 2000, Syphard et al. 2011).The Biomass Harvest extension selects and removes biomass based on specific timber harvest prescriptions.Stands are selected for harvest based on a user defined ranking (e.g., when cohorts reach economically merchantable age).Within a selected stand, user defined prescription criteria are applied, including patch size, percent biomass removed, whether harvests are repeated, species specific removal, planting, etc. Rotation periods define the amount of the landscape harvested within each prescription and management area for each time step (Gustafson et al. 2000).
To simulate BAU management, we simulated current harvesting regimes practiced in northeastern Minnesota and northern lower Michigan.We delineated stands within management areas based on ownership groups (Fig. 1).In the Minnesota landscape, we utilized landowner group specific BAU harvesting regimes previously described (Ravenscroft et al. 2010).We adjusted these prescriptions based on expert opinion (A.W. D'Amato, personal communication) (see Appendix: Table A1).In the Michigan landscape, we included specific prescriptions used by Michigan DNR and the Huron-Manistee National Forest and made adjustments for Private Non-Industrial Forest Land (Potter-Witter 2005) and Private Industrial Land (C.Webster, personal communication) (see Appendix: Table A2).Our BAU harvesting regimes represent our best guess at currently practiced management.They do not portray precisely what is nor will be practiced in the future.
We defined management areas (i.e., ownership groups) using both the Conservation Biology Institute (CBI) Protected Areas Map (The Conservation Biology Institute 2010), and the FIA program forest ownership map (Nelson et al. 2010).We used a combination of the two maps as the CBI map did not delineate private industrial land from private non-industrial land but otherwise had better land owner precision.Our management area delineation assumes that ownerships will stay static within aggregated groups.Within management areas, we implemented the individual prescription rotation periods based on the proportion of each prescription cover type present at year zero.We delineated individual stands within management areas based on the IFMAP forest types in Michigan and used previously delineated stands from northeastern Minnesota (Ravenscroft et al. 2010).
Along with individual species and species richness, we considered species assemblages as representative of forest types (Eyre and Society of American Foresters 1980, Bøhn and Amundsen 2004) (see Appendix: Table A3).An analysis of forest type provides summary information about species groups.To quantify spatially explicit species Aboveground Biomass (AGB) (g m À2 ) and to classify forest type, we utilized the LANDIS-II Biomass Output extension (version 2.0) and the Biomass Reclassification extension (version 2.0) respectively.Although novel forest types are likely to emerge in the face of climate change (Williams et al. 2007), we focused on the current suite of forest types in our assessment of the effects of climate change.We identified forest types using indicator species common on the landscape but relatively unique to each forest type.We adjusted forest reclassification species assignment based on matching the proportion of individual forest types found in regional FIA plots to the proportion of individual forest types found in LANDIS-II cells at year zero (see Appendix: Table A4).To process outputs and visualize results, we used the raster library in R (R Development Core Team 2011).

Data analysis
We initially used five replicates to explore the stochastic variation within scenarios caused by the stochastic nature of disturbance events.Maximum variance of AGB at year 2150 was ,2% of the mean for each landscape.Given the low among-scenario variance caused by the low stochastic variation in disturbance events compared to the size of the landscapes, we selected one replicate from each climate scenario and landscape for further analysis.For each tree species, we calculated AGB as the mean biomass density (g m À2 ) across the landscape for each scenario and time step.We calculated AGB total on each landscape as the total above ground biomass (Tg) on each landscape for each scenario and time step.We considered AGB total in order to assess each landscapes biomass contribution to the region.
To assess species diversity, we calculated Shannon index of diversity (H 0 ) using AGB abundance of tree species.In order to reduce the magnitude of inaccuracies comparing diversity across scenarios, we transformed Shannon indices to the effective numbers of species (e H 0 ) (Jost 2006).Effective number of species is the number of species present if all species were equal in abundance.We calculated diversity within individual cells across each landscape for each climate change scenario at 10-year time steps starting after year 2050.As a measure of productivity within each ecoregion, we used aboveground annual net primary productivity (ANPP) as calculated within the Biomass Succession extension.
As our simulations did not account for migration of new species from outside the landscapes, we assessed diversity within a 10km buffer from the southern edge of each landscape.Scheller and Mladenoff (2008) found movement of common southern Great Lake species less than 10 km within 200 years.This rate is comparable to molecular studies of species migration under climate change (McLachlan et al. 2005).The Minnesota landscape is largely surrounded by boreal lowland forest and Lake Superior.The Michigan landscape is largely surrounded by Lake Michigan, and Lake Huron (Fig. 1).Additional limits to species migration into our landscapes include fragmentation of forests due to agriculture, and competitive lags to establishment.
For each simulation, we assessed the relationship between diversity and productivity at ecoregion levels using linear regression models.To assess landscape variability of diversity, we calculated average cell diversity within ranked climate and soil regions.Using management area maps (Fig. 1), we also calculated mean diversity within ownership groups.All data analysis was completed using R (R Development Core Team 2011).Full model validation is not possible against future empirical data.Nonetheless, to evaluate our ability to simulate AGB for initial conditions, we analyzed AGB estimates as measured in the 2140 FIA field plots compared to the v www.esajournals.orgimputed LANDIS-II year zero biomass using the Pearson's correlation coefficient and the root mean squared error (RMSE).While this method does not truly validate future successional trajectories, a strong relationship between empirical estimates and modeled biomass after the spin-up phase suggests that the spin-up is comparable to current successional trajectories of biomass (Thompson et al. 2011).In addition, most of the species-specific parameters, PnET-II, LANDIS-II, and the Biomass Succession extension have been evaluated in similar landscapes in the Great Lakes Region (Scheller and Mladenoff 2004, Radeloff et al. 2006, Xu et al. 2007, Ravenscroft et al. 2010).

RESULTS
For each FIA plot imputed onto the northern lower Michigan landscape, our evaluation of AGB resulted in a Pearson's correlation between initial (year 2010) LANDIS-II AGB and observed FIA estimated AGB of 0.64.The RMSE, representing the spread of the fit, was 44.3 Mg ha À1 .

Productivity and diversity
Across both Minnesota and Michigan landscapes, we found weak to no diversity-productivity relationships under current and low emission climate scenarios (Fig. 5).Under these scenarios, effects of soil water holding capacity (SWHC) stratify ecoregions along a productivity v www.esajournals.orggradient.Low SWHC was associated with lower productivity.Average simulated productivity across ecoregions at year 2150 was similar for current climate (206 and 236 g m À2 yr À1 ) and the low emissions scenario (216 and 227 g m À2 yr À1 ) for the Minnesota and Michigan landscapes, respectively.Under the high emissions climate scenario, however, we found positive relationships between diversity and productivity in both landscapes.Linear regression slopes were 11.9 and 21.1; adjusted r 2 values were 0.043 and 0.398 for the Minnesota and Michigan landscapes, respectively.Within this climate scenario, average simulated productivity at year 2150 was 87 and 82 g m À2 yr À1 , representing 58% and 65% lower ANPP compared to current climate in the Minnesota and Michigan landscapes, respectively.

Landscape variation in diversity
Simulations in the Michigan landscape resulted in higher tree species diversity compared to the Minnesota landscape (Fig. 6).Under all climate scenarios, high SWHC did not result in the highest simulated diversity within each landscape, however the lowest SWHC regions resulted in the lowest diversity under the climate change scenarios.Within the Minnesota land- v www.esajournals.orgscape, the warmest climate region resulted in the highest diversity under the current climate scenario.Under the high emissions climate scenario, the coolest climate region had the highest diversity.In the Michigan landscape, the warmest climate region resulted in the lowest diversity under both climate change scenarios.In the Minnesota landscape, the forest reserve management areas had the highest diversity for all three climate scenarios and all time steps.In the Michigan landscape, private industrial forests (PIF) had seemingly high diversity under the current and low emission climate scenario.These lands represent a very small proportion of the landscape and with a small number of cells could easily have resulted in low diversity given the management history and relative size of the landscape (Fig. 1).In the Minnesota landscape, diversity across climate scenarios resulted in the highest diversity in the low emissions climate scenario.In the Michigan landscape, diversity across climate scenarios resulted in higher diversity in the current climate, and progressively lower diversity in the low and high emissions climate scenario.

Aboveground biomass
The initial simulated AGB total in the Michigan landscape was substantially larger than the Minnesota landscape (Fig. 7) due in large part to the larger Michigan landscape (Fig. 1).Simulated forest growth resulted in an initial increase for all climate scenarios in both landscapes.For both landscapes, the AGB total in the low emissions climate scenario resulted in a trajectory almost identical to that of current climate.AGB total in the high emissions climate scenario resulted in a departure after approximately two decades and decreased to approximately initial levels by year 2150.
We found substantial differences in AGB for individual species among landscapes and their response to climate change.Most species had higher AGB in the Michigan compared to the Minnesota landscape (Fig. 8), except for black spruce (Picea mariana), quaking aspen, balsam fir, and paper birch.Under current climate, the simulated species response suggests that both landscapes would favor more late successional species such as balsam fir and white spruce (Picea glauca) in the future.The response to climate change suggests that boreal species at the southern edge of their range will decline in both Minnesota and Michigan landscapes.Compared to current climate, both climate change scenarios generally resulted in less black spruce, white spruce and paper birch (Betula papyrifera) AGB in both landscapes.In contrast, many northern hardwoods, including red maple, sugar maple (Acer saccharum) and American basswood (Tilia americana), increased rapidly in response to climate change for the first 80 years of our simulations.

Forest type shifts
In the Minnesota landscape, the current climate scenario under BAU management resulted in an increase in the late successional spruce-fir   A1FI) assuming BAU forest management.American beech, black oak, northern pin oak, white ash, balsam poplar, eastern hemlock and Scotch pine were not included in the Minnesota landscape and bur oak was not included in the Michigan landscape due to low or no abundance.Order of species (left to right, top to bottom) is based on rank species biomass abundance at year zero in the northeastern Minnesota landscape.v www.esajournals.orgforest type between 2000 and 2150 (Fig. 9).The aspen-birch and red/jack pine forest types slightly declined.The northern hardwoods, white pine, and oak forest types all persisted through the simulation as small components of the landscape.The climate change scenarios resulted in a departure from current conditions as northern hardwoods and white pine forest types became more dominant on the landscape and the sprucefir type declined.The climate change scenarios resulted in a larger magnitude of change under the high emission compared to the low emission scenario.Under both climate change scenarios, the Boundary Waters Canoe Area Wilderness provided a refuge for the spruce-fir forest type (Fig. 9; northern component of the Minnesota landscape).
In the Michigan landscape, the current climate scenario under BAU management resulted in substantial increases in the white pine forest type and minor shifts for other forest types (Fig. 10).Both climate change scenarios resulted in an increase in northern hardwoods and oak association forest types and a decrease in spruce-fir.By simulation year 2150, the low emission scenario resulted in a greater increase in white pine and a smaller increase in oak forest types compared to the high emissions scenario.
Our results suggest that climate change may alter the relationship between diversity and productivity under conditions of declining productivity.Furthermore, diversity may be affected by regional climate differences, soils, and forest management.In addition, our results suggest less total biomass under a high emissions climate scenario compared to current climate.

DISCUSSION
With the threat of climate change, there is v www.esajournals.orgstrong support to preserve and manage for biodiversity in the ecological literature (e.g., Chapin et al. 1998, Hampe and Petit 2005, Wilson 2010) as well as the forest management literature (e.g., D'Amato et al. 2011, Puettmann 2011).A future decline in biodiversity has the potential to decrease productivity (Tilman et al. 2012).The positive diversity-productivity relationship under the high emission climate scenario, and the lack of relationship under current and low emissions climate scenarios (Fig. 5) supports the diversity-productivity hypothesis that diversity will follow a transposed bell shaped curve with productivity (Loreau et al. 2001, Whittaker et al. 2001).Under this hypothesis, a positive diversityproductivity relationship exists under relatively low productivity, and a negative diversity-productivity relationship exists under relatively high productivity.As simulated productivity decreased with climate change, competitive exclu-sion (the ability of one species to dominate over another) decreased as complementarity (the ability of a diverse community to use resources more completely) increased.Regardless of the cause and effect relationship, the association between diversity and productivity may become stronger under less favorable conditions due to climate change.Whittaker et al. (2001) describe at what spatial and temporal scale productivity and diversity may be more regulated by top-down explanatory variables compared to bottom-up local influences.Our results imply a bottom-up regulated system under current and low emissions climate where productivity and diversity is largely regulated by local variables (e.g., SWHC).Under a high emissions climate future, however our results indicate more of a top-down regulated system.Climate change may impinge on productivity-diversity relationships at a finer spatial resolution than under current climate.

Contrasting landscapes
Our results indicate that the Michigan landscape has higher diversity of tree species (Fig. 6) and forest types (Fig. 10) compared to the Minnesota landscape.These results reflect biotic, abiotic, and management differences.The pattern of diversity within forest reserves highlights the difference in overall management in the Minnesota and Michigan landscapes (Fig. 6).The reserves in the Minnesota landscape result in a larger contrast to the actively managed forests compared to the Michigan landscape.Within the actively managed forests, the Minnesota landscape includes more frequent low diversity aspen silviculture systems (D'Amato et al. 2009) compared to the Michigan landscape (see Appendix: Table A1 and A2).In the Minnesota landscape, 50% of the landscape was initially classified as aspen-birch (Fig. 9) and reflects the large component of aspen silviculture in the managed landscape.In the Michigan landscape, forest types were more equally represented.Aspenbirch accounted for only 22% of the landscape (Fig. 10).
These results suggest that the diversity of forest types and species in the Michigan landscape may be better poised to adapt to climate change compared to the Minnesota landscape.ANPP declines more in the Michigan landscape under high emissions climate than the Minnesota landscape, however.Although more diverse than Minnesota, the suite of species simulated under BAU management in Michigan and the high emissions climate is not capable of maintaining the equivalent high productivity simulated under current climate.Our simulations, however, did not allow the migration of novel species from outside the landscape.Future forest resilience may depend on the ability of novel species to migrate at the pace of changing climate (Buma and Wessman 2013).
In addition to higher diversity in less intensively managed forests, our results indicate a general pattern of higher diversity within cooler climate regions and more productive soils (Fig. 6).Regional variation in climate may have profound effects on species composition in the Great Lakes region (Davis et al. 2000) and beyond (Davidson andJanssens 2006, Svenning andSandel 2013).For example, Davis et al. (2000) found, in 10,000 years of pollen records, that species composition changed less in sites within close proximity to Great Lakes where temperature and precipitation were moderated.In the Minnesota landscape, diversity declined along a gradient of progressively warmer climate regions under the high emissions climate scenario.However, the Michigan landscape resulted in a less variation in diversity due to climate region.This may be due to the inherent higher diversity within all climate regions in the Michigan landscape.
Notwithstanding climate, available soil water and soil type in general, may affect future species diversity (Knapp et al. 2002, Anderson andFerree 2010).In the Michigan landscape, there was a positive relationship between SWHC and diversity in all climate scenarios.In the Minnesota landscape, the effect of SWHC on diversity was less distinct.Although we did not explore this hypothesis, regions of high SWHC may be currently managed for intensive even-aged aspen silviculture which would result in lower diversity.

Aboveground biomass
There is increasing interest in forests contribution to carbon sequestration and storage (Lal 2005, Fahey et al. 2010, Bo ¨ttcher et al. 2012, Clemmensen et al. 2013).Under all climate futures, our results indicate an initial large increase in AGB total .The capacity for increased AGB in both landscapes is likely due to forest recovery following a history of intense timber harvesting (Pyne 1982, Flader 1983, Bo ¨ttcher et al. 2012).As forests mature, the sustained increase in AGB total expected under current climate can be further explained by the projected increase in canopy structural complexity (Hardiman et al. 2011).In addition, the simplification of the initial communities within the modeling framework results in a simulated increase in cohort heterogeneity.Under the high emissions climate scenario our simulations indicate a decline in AGB (Fig. 7).Empirically, Fisichelli et al. (2012) studied saplings across a boreal-temperate transition gradient in Minnesota and found that with increasing temperature, height and radial growth rates of boreal species declined and those of temperate species increased.
Successful insurance species are capable of responding positively to changing conditions (Naeem andLi 1997, Walker et al. 1999).Our results suggest that the species expected to be insurance species (i.e., northern hardwoods) especially in the warmer southern landscape, may not be tolerant of a future warmer climate projected under the high emissions scenario.Poorly suited species under a high emissions climate future has been documented in other regions (Buma and Wessman 2013).Our simulations result in declining AGB of black spruce, white spruce, balsam fir, jack pine, and paper birch under climate change.Common trends of increasing species abundance under climate change include white pine (Pinus strobus), black cherry (Prunus serotina), American basswood (Tilia americana), and American elm (Ulmus americana).All four of these species from the northern lower Michigan landscape, however, indicate a decline before the end of the simulation (Fig. 8).The species that share a simulated increase or decrease across landscapes are consistent with empirical observations (Fisichelli et al. 2012, Fisichelli et al. 2013) and other modeling throughout the Great Lakes region (Iverson et al. 2008, Xu et al. 2009, Ravenscroft et al. 2010).
A climate change increase in AGB total composed of responding hardwood and oak species was found by other modeling research in the region (Scheller and Mladenoff 2005).In contrast, both of our landscapes had declining AGB total under the high emissions scenario after year 2050 compared to the current climate scenario (Fig. 7).Compared to current climate, the high emissions climate scenario resulted in an initial increase in the hardwood and oak species (e.g., red maple, American basswood, black cherry, northern red oak [Quercus rubra]).This increase, however, did not make up the difference in lost AGB by boreal species.The high emissions climate future that we used (i.e., GFDL A1FI) includes higher CO 2 and temperature projections than the A2 climate futures used by Ravenscroft et al. (2010) or the IS92A climate future used by Scheller and Mladenoff (2005) (IPCC 2007).In addition, our work simulated an increasing trend of climate warming beyond year 2100 while previous work held climate constant after year 2100.

Forest type
Ecotones around the globe are valuable and vulnerable to the effects of climate change (Smith et al. 1997, Goldblum andRigg 2010).The Minnesota landscape is centered closer to the boreal-temperate ecotone while the Michigan landscape is located further into the temperate forest matrix (Albert 1995).Our results suggest that the Minnesota landscape will experience rapid forest type conversion under the high emissions climate scenario (Fig. 9) as expected within an ecotone (Allen and Breshears 1998).These results corroborate previous forest simulation research indicating forest type conversions expected under climate change in the region (Iverson and Prasad 2001).The Michigan landscape also indicates rapid forest type conversion, however the large forest type shifts occur under all three climate scenarios (Fig. 10), not only within the low and high emissions climate futures.This suggests that forests in the Michigan landscape (more than the Minnesota landscape) have strong successional trajectories expected, notwithstanding climate changes.Our simulated transition of early to more mid-successional species in both landscapes under current climate is corroborated by observed shifts in northern lower Michigan (Gough et al. 2010) and Minnesota (Hanberry et al. 2012).
The future of early successional forests will depend on highly variable disturbance regimes, which are difficult to predict.The homogenizations of aspen forests are closely linked to harvesting; our results reflect future harvest regimes based on our best approximation of BAU management.Harvesting regimes are largely a product of market demand for different wood products.The recent recession has resulted in a decline in aspen harvests for pulpwood (Woodall et al. 2011).If this trend continues, we expect to see a shift to more shade tolerant species.A change in future markets however could change harvesting regimes and resulting forests, considerably.Further uncertainties may also include changing regulations, technology, and landowner priorities.
Forest type changes in the Minnesota landscape indicate more of a climate effect compared to the Michigan landscape.Hooper and Vitousek (1997) experimentally found that both species composition and the number of functional groups positively affected grassland productivity.The higher diversity of species biomass (Fig. 6) and of forest types in the Michigan compared v www.esajournals.orgto Minnesota landscapes may contribute to the greater adaptive capacity of future forest types.

Climate and management uncertainty
Our projections assumed that climate warming will continue beyond year 2100 based on trends and variability projected by GCM's (IPCC 2007) and downscaled climate (Stoner et al. 2012).In addition to the large uncertainty in both the GCM and emission scenarios, our results should be interpreted with a gradient of increasing uncertainty beyond the range of GCM projections.Although observed carbon emissions suggests that our current carbon emission trajectory is at or above that of the high emissions scenario (Raupach et al. 2007, Jennings 2012, Peters et al. 2012), the low emissions scenario coupled in our simulations to a less sensitive GCM indicates substantial differences in ANPP, AGB, and diversity from the high emissions scenario.Nevertheless, even with current climate futures, successional trajectories indicate the potential for future novel ecosystems rather than a stable state.
Although our analysis of diversity a 10-km buffer at the southern edge of each landscape to account for potential species migrating north, species migration rates vary and could exceed our 10-km spatial restriction (Scheller and Mladenoff 2008).Our simulations assumed constant disturbances (harvesting, wind, and fire) over time.More frequent disturbances could allow more adaptation potential (Buma and Wessman 2013).Nevertheless, forest distribution shifts are expected to be slow and will be limited by fragmentation, competition, and the island-like design of our landscapes (McLachlan et al. 2005, Vanderwel andPurves 2013).
Our results should not be interpreted as predictions.Rather, scenario projections illuminate interactions between processes and regional trends.And like all simulation models, our results represent a simplification of reality.There are processes that we did not include in the modeling framework but that likely will be important.These include, but are not limited to the effects of migrating invasive earth worms which may change the soil substrate required for regeneration (Holdsworth et al. 2007, Hale et al. 2008, Larson et al. 2010) and the effects of browse damage by large ungulates (Anderson et al. 2002, Danell et al. 2003, Frelich and Reich 2010, Fisichelli et al. 2012).If deer browsing continues to affect forest communities at present rates, the affects will be severe and extensive.
The effects of CO 2 fertilization, without considering additional limits to growth, will increase the photosynthetic rate of future forests (Mickler et al. 2002, Norby et al. 2005).Growth enhancement by CO 2 fertilization will be limited by Nitrogen (Ko ¨rner 2006, Hyvo ¨nen et al. 2007, Reich and Hobbie 2012), which we did not include as a dynamic process.Additional research suggests that compared to the forest response to prior land use history, CO 2 fertilization has little affect (Caspersen et al. 2000).Furthermore, we did not include the influence of ozone, which will likely diminish the effects of CO 2 fertilization (Wittig et al. 2008, Ainsworth et al. 2012) and may decrease overall growth especially in northern lower Michigan where ozone pollution is more prevalent (Koo et al. 2012).Globally, no increase in NPP has been observed by satellite measurements for the last 30 years of increased atmospheric carbon dioxide (Running 2012).
In addition to emission and climate uncertainty, there exists tremendous uncertainty about how future forests will be managed.Our simulations assume BAU timber harvest which includes the continuation of intensive resource extraction.Current timber management relies heavily on the continued production of fiber for the pulp industry (Hanberry et al. 2012).Additional modeling research should address the growing trend (D'Amato et al. 2009) and interest (Millar et al. 2007) to manage forests for complexity (Puettmann 2011) and adaptability (Kuparinen et al. 2010).In addition, there is support to preserve the trailing edge of a migrating ecotone which may contain high genetic diversity (Hampe andPetit 2005, Willis et al. 2010).Future research will address management's ability to affect the resistance and resilience of future forests to the effects of climate change.
Megan Creutzburg, David Neumann, Klaus Puetmann, and Christopher Webster for their generous assistance, review and advice.Additionally, we thank two anonymous reviewers for insightful comments and perspectives.This project was funded by the U.S. Fish & Wildlife Service Upper Midwest and Great Lakes Landscape Conservation Cooperative.Additional funding was provided by the Portland State University Department of Environmental Science and Management.

Fig. 1 .
Fig. 1.Inserted map of study landscapes within continental United States.Ownership groups within study landscapes in northeastern Minnesota and northern lower Michigan.

Fig. 2 .
Fig. 2. Nested climate and soil regions used to create ecoregions in northeastern Minnesota (left) and northern lower Michigan (right).Climate regions are designated by numbered polygons in each landscape.Nested soil regions represent binned SWHC (cm/m depth) values and are designated by colored pixels within climate regions.Note that finer resolution of SWHC data were available in northern lower Michigan and not northeastern Minnesota.

Fig. 4 .
Fig. 4. Aggregated climate regions based on mean July temperature in each landscape.Labels indicate average current July temperature (8C) from 30 years of PRISM observations.Note that warmest climate region in Minnesota landscape is cooler than coolest climate region in Michigan landscape.

Fig. 5 .
Fig. 5. Relationship between diversity and aboveground ANPP (g m À2 yr À1 ) at year 2150 in northeastern Minnesota (top) and northern lower Michigan (bottom) for three climate scenarios.Unique point symbols represent soil regions across a gradient of SWHC values (cm/m).Diversity is represented by effective number of species (e H 0 ) and is calculated from Shannon index.

Fig. 6 .
Fig. 6.Diversity (e H 0 ) across the last 100 years of simulation (i.e., 2050 to 2150) in northeastern Minnesota and northern lower Michigan.Each group of clustered data points represents a single time series plot for a given climate scenario.Variation in diversity within soil water holding capacity (SWHC) regions, climate regions, and management areas is represented by different colored and shaped points.Climate region plots represent diversity within aggregated climate regions (climate region 1 ¼ warmest climate region, climate region 5 ¼ coolest climate region).Management area plots represent diversity within forest management groups.DNR represents both DNR and county lands in Minnesota and only DNR lands in Michigan.

Fig. 7 .
Fig. 7. Total aboveground biomass (AGB total ) within the northeastern Minnesota and northern lower Michigan landscapes under three climate scenarios assuming BAU management.

Fig. 8 .
Fig. 8. Aboveground biomass (AGB) (g m À2 ) by species for 150 simulation years across the northeastern Minnesota (MN) and northern lower Michigan (MI) landscapes for three climate scenarios (Current, B1, and A1FI) assuming BAU forest management.American beech, black oak, northern pin oak, white ash, balsam poplar, eastern hemlock and Scotch pine were not included in the Minnesota landscape and bur oak was not included in the Michigan landscape due to low or no abundance.Order of species (left to right, top to bottom) is based on rank species biomass abundance at year zero in the northeastern Minnesota landscape.

Fig. 9 .
Fig. 9. Northeastern Minnesota map of classified forests type at simulation year 2000 (upper left), proportion of forest types changing through time for each climate future (middle), and associated forest type maps at simulation year 2150 (right).Figure colors consistent throughout.

Fig. 10 .
Fig. 10.Northern lower Michigan map of classified forest type at simulation year 2000 (upper left), proportion of forest types changing through time for each climate future (middle), and associated forest type maps at simulation year 2150 (right).Figure colors consistent throughout.
. Climate change projections for the GCMs and emission scenarios we utilized are summarized in Handler et al. (in press a) and Handler et al. (in press b) for Minnesota and Michigan, respectively.For each unique climate region, we accessed downscaled temperature and precipitation data through the USGS data portal http://cida.usgs.gov/climate/gdp/(Stoner et al. 2012).To estimate projected PAR, we accessed point source projected radiation data from meteorological stations throughout the landscape.After observing the low variability in PAR across observation sites, we used mean monthly PAR within each climate region.