Introduction

The Aspen-FACE experiment was an 11-year study of the effect of elevated CO2 and ozone (alone and in combination) on the growth of model aspen communities (pure aspen, aspen-birch, and aspen-maple) in the field in northern Wisconsin, USA. Uncertainty remains about how these short-term plot-level responses might play out over broader temporal and spatial scales where climate change, competition, succession, and disturbances interact with tree-level responses. In this study, we used a new physiology-based approach (PnET-Succession v3.1) within the forest landscape model LANDIS-II to extrapolate the FACE results to broader temporal scales (and ultimately to landscape scale) by mechanistically accounting for the globally changing drivers of temperature, precipitation, CO2, and ozone. We added novel algorithms to the model to mechanistically simulate the effects of ozone on photosynthesis through ozone-induced impairment of stomatal control (i.e., stomatal sluggishness) and damage of photosynthetic capacity at the chloroplast level.

Results

We calibrated the model to empirical observations of competitive interactions on the elevated CO2 and O3 plots of the Aspen-FACE experiment and successfully validated it on the combined factor plots. We used the validated model to extend the Aspen-FACE experiment for 80 years. When only aspen clones competed, we found that clone 271 always dominated, although the ozone-tolerant clone was co-dominant when ozone was present. Under all treatments, when aspen clone 216 and birch competed, birch was always dominant or co-dominant, and when clone 216 and maple competed, clone 216 was dominant, although maple was able to grow steadily because of its shade tolerance. We also predicted long-term competitive outcomes for novel assemblages of taxa under each treatment and discovered that future composition and dominant taxa depend on treatment, and that short-term trends do not always persist in the long term.

Conclusions

We identified the strengths and weaknesses of PnET-Succession v3.1 and conclude that it can generate potentially robust predictions of the effects of elevated CO2 and ozone at landscape scales because of its mechanistically motivated algorithms. These capabilities can be used to project forest dynamics under anticipated future conditions that have no historical analog with which to parameterize less mechanistic models.

Global changes to atmospheric composition have the potential to impact forest productivity, succession, and resilience. Carbon dioxide concentration (hereafter denoted as simply CO2) has increased steadily since the beginning of the industrial revolution because of land clearing and fossil fuel burning and it is widely expected to continue increasing well into the current century (Siegenthaler et al. 2005). Combustion of fossil fuel also releases nitrogen oxides (NOx), which react with O2 and volatile organic compounds, primarily methane (CH4), in the presence of sunlight to produce ozone (O3). Increased CO2 is known to increase plant productivity and has stimulated growth of aspen forests in the upper Great Lakes region of the USA over the last 50 years (Cole et al. 2010). Such changes in productivity may not affect all species or all canopy positions equally, which may be altering composition of some terrestrial ecosystems (Bond and Midgley 2000; Kgope et al. 2010). On the other hand, tropospheric O3 is considered the most significant air pollutant negatively affecting forest productivity worldwide (Matyssek et al. 2013). Ozone is a potent oxidizer that disrupts cell function, resulting in decreased plant productivity, leaf necrosis, and increased mortality (Karnosky et al. 2005). Although tropospheric ozone concentration (hereafter denoted as O3) has generally decreased in the continental USA since 2002, some regions in the USA and globally have experienced increases in O3 production because of increased emissions of chemical precursors (Cooper et al. 2012). Thus, a number of forested regions continue to be subjected to phytotoxic levels of O3.

Many experiments on the effects of elevated CO2 and O3 on tree growth have been conducted at growth chamber and field plot scales (Kozovits et al. 2005; Hoshika et al. 2012a). Free-Air Carbon Enrichment (FACE) experiments represent the most robust attempts to demonstrate the effects of elevated CO2 on trees in a realistic field setting (Karnosky et al. 2001; Paoletti et al. 2017). The Aspen-FACE experiment was unique in that it experimentally modified both CO2 and O3 concentrations in forest ecosystems in a replicated, factorial experiment, generating a record of the effects of those gases on the growth, health, and mortality of native aspen (Populus tremuloides), paper birch (Betula papyrifera), and sugar maple (Acer saccharum) in northern Wisconsin, USA (Kubiske et al. 2007). These tree species represent some of the most economically and ecologically important species in northern hardwood and sub-boreal forests of eastern North America (Dickson et al. 2000), and the management of the landscapes in which they are found presents critical challenges in the face of rising pollution from fossil fuel combustion (Smith 2012). However, results from such short-term, local studies must be extrapolated to longer and broader scales to usefully inform the management of ecosystems, and techniques to conduct this scaling are rudimentary at best.

In most modeling studies (e.g., Hanson et al. 2005), the effects of CO2 and O3 on photosynthesis are separately modeled and do not account for the interactive effect of CO2 and O3. In addition to the Aspen-FACE experiment, several studies found that elevated CO2 ameliorated O3 reductions to photosynthesis (Grams et al. 1999; Kitao et al. 2015). In addition to causing tissue damage, ozone reduces stomatal function, generally resulting in an increase of CO2 in leaf interiors, a slower stomatal response to environmental stimuli (Paoletti and Grulke 2005) and increased water loss. Modeling schemes are needed that account for these O3 effects on photosynthesis, including the interactions with elevated CO2, especially when simulating future scenarios where environmental conditions have no historical analog (Ito 2010).

There has been one attempt to scale Aspen-FACE results to the stand scale (Karnosky et al. 2005), and one attempt to scale them to the landscape scale (Gustafson et al. 2013), but uncertainty remains high about how these short-term plot-level responses might play out at the landscape scale when climate change, elevated CO2 and O3, competition, succession, and disturbances all interact with tree-level responses over broader temporal and spatial scales (i.e., centuries and landscapes, respectively). For example, Gustafson et al. (2013) used Aspen-FACE data to calibrate the Biomass Succession extension of the forest landscape model LANDIS-II (Scheller et al. 2007) to the results of the FACE experiment and project the experiment into the future at both site and landscape scales. However, growth processes in the Biomass Succession extension are simulated using a relatively phenomenological approach, where maximum aboveground net primary productivity (ANPP) of a species is estimated as an average under specific climate and atmospheric conditions. Although the parameter Maximum-ANPP can be varied through time, it does not respond to independent variation in specific drivers (e.g., precipitation, temperature, CO2, O3) at a sub-decadal temporal scale, making it insensitive to extreme events. It is becoming clear that extreme events such as long droughts can have a major effect on forest dynamics (Clark et al. 2016), and not simulating these events reduces confidence in the projections. Furthermore, Biomass Succession does not explicitly include processes of competition for light and water, instead simulating competition for “growing space,” assuming adequate light and water for all cohorts of all species on a site. Most importantly, because Biomass Succession uses a phenomenological approach, its parameterization for novel conditions that have never been empirically studied introduces considerable uncertainty (Gustafson 2013).

In this study, we used a more mechanistic LANDIS-II succession extension (PnET-Succession v3.1) that uses a “big leaf,” physiology-based approach based on first principles to scale the Aspen-FACE results to landscape scales by mechanistically accounting for the effects of CO2 and O3 on photosynthesis and competition. Our objectives were (1) to generate a more scientifically plausible simulation of the Aspen-FACE experiment (i.e., calibrate and validate a more mechanistic model) and (2) to apply the validated model to extrapolate the empirical results to conditions not covered by the empirical data (i.e., longer time periods and other assemblages of tree species). Our ultimate objective is use the LANDIS-II model to project the long-term effects of elevated CO2 and O3 on the dynamics of northern hardwood and sub-boreal forests of North America (including seed dispersal and disturbance), which will be completed in a subsequent study. This study is designed to function as a “proof of concept” to assess whether the more mechanistic approach can mimic empirical observations and increase confidence in model predictions. Such mechanistically derived results can reduce uncertainty about how future atmospheric and climate conditions might impact forest dynamics and the future abundance of species that have economic and ecological value, particularly in regions prone to harmful levels of ozone.

Our study is built on a few major assumptions. First, our purpose was to build a modeling capability to robustly project elevated CO2 and O3 effects to landscape spatial and temporal scales, able to include other disturbances and factors (e.g., seed dispersal, climate) that determine forest ecosystem dynamics. Because future CO2 and O3 concentrations are expected to range outside of historically observed levels, a mechanistic modeling approach based on physiological first principles is preferable to phenomenological approach based on behavior observed under past conditions (Gustafson 2013). Because our model must be tractable at landscape spatial and temporal scales, we must use some phenomenological components in an otherwise mechanistic model, which is likely to produce some discrepancies with plot-level observations (Aspen-FACE experiment), but the magnitude of the discrepancies is expected to be small. Secondly, we assumed that the Aspen-FACE experiment can be generalized to other places and other species, which was a key objective of that experiment (Dickson et al. 2000). Our study is an attempt to generalize these experimental results. Finally, PnET-Succession assumes that cohort growth and competition is primarily driven by foliar nitrogen (surrogate for innate photosynthetic capacity (Amax)), modified by the limiting factors of water, light, and temperature. In our study, we also assumed that concentrations of CO2 greater than 365 ppm partly enhance photosynthetic capacity by reducing stomatal conductance (as in PnET-II), which reduces water loss and ozone absorption. Ozone is assumed to be a similar limiting factor, decreasing Amax. In addition, it is assumed that ozone impairs the ability of stomata to control balance of carbon gain and water loss (i.e., stomatal sluggishness (Paoletti and Grulke 2005)). Competition can differentially alter light and water availability among cohorts on a site, but competition does not directly alter exposure to CO2 and O3 in our model. Thus, direct effects of CO2 and O3 on competitive interactions between cohorts in the real world are not included in the model, but a number of indirect effects can readily occur. All these assumptions were held constant across treatments, so differences in responses can be attributed to the treatments, given the initial conditions and the assumptions.

Model description

LANDIS-II is a process-based forest landscape disturbance and succession model that simulates the forest development processes of dispersal, establishment, growth and competition, and the forest degenerative processes of senescence and disturbance at large spatial (> 105 ha) and temporal (centuries) scales. Landscapes are represented as a grid of interacting cells with user-defined cell size. Individual cells are assumed to have homogeneous canopy layering and are spatially aggregated into land types with homogeneous climate and soils. Forest composition at the cell level is represented as age cohorts of individual tree species that interact via a suite of vital attributes (e.g., shade tolerance, fire tolerance, seed dispersal, ability to sprout vegetatively, longevity) to produce nondeterministic successional pathways sensitive to disturbance type and severity (Mladenoff 2004). We used LANDIS-II version 6.2 (Scheller et al. 2007), which consists of a core collection of libraries and a collection of extensions that represent specific ecological processes.

To simulate the growth and competition processes, we used the PnET-Succession extension (De Bruijn et al. 2014) that uses algorithms from the PnET-II ecophysiology model (Aber et al. 1995) to simulate cohort growth and competition as a function of competition for light and water. PnET-Succession upscales small scale (i.e. ,per gram foliage) biochemical processes such as photosynthesis, respiration, and transpiration to the grid-cell using a “big leaf” approach by integrating light extinction and water consumption in stacked canopy and sub-canopy layers and computing a dynamic soil water balance. First, species-cohort photosynthetic capacity under optimal conditions (Amax) is determined as a function of foliar nitrogen concentration. Actual photosynthesis (A) in a given month and stratum of the canopy is computed by applying reduction multipliers (0.0–1.0) that reflect departure from optimal conditions (stress) for several factors. The primary factors are light and soil water availability, but temperature, age (senescence), and ozone dose are also included. Soil water is tracked at the grid-cell level using a simple bulk-hydrology model based on precipitation, loss to evaporation, runoff and percolation out of the rooting zone, and transpiration by species cohorts. Cohorts compete for water and light in each cell, and access to light and soil moisture is proportional to cohort biomass, with some stochasticity introduced to give smaller cohorts some access to these resources when they are limited. When water is adequate, the rate of photosynthesis for a given species cohort is determined by the light that is available to the cohort (dependent on canopy position and leaf area), and decreases with age and departure from optimal temperature. As soil water availability decreases, photosynthesis also decreases according to the water stress reduction factor (fWater), such that fWater = 0.0 when water potential is less than the species’ drought tolerance (H4, Table 1). PnET-Succession accounts for reductions in photosynthesis due to growth and maintenance respiration using a Q10 relationship (Atkins 1978), such that foliar respiration rate depends on temperature and moisture, while maintenance respiration depends on temperature and biomass. Net photosynthetic production is allocated to biomass pools of foliage, wood, root, and reserves (non-structural carbon) according to allocation parameters (FracFol, FracBelowG, dNSC; Tables 1 and 2). Establishment of new cohorts is stochastic with species-specific establishment probabilities calculated monthly as a function of soil water and sub-canopy light. Other processes that kill cohorts or remove some of their biomass are simulated by independent disturbance extensions. The interaction of grid cells via dispersal and spatial disturbance processes within the LANDIS-II framework robustly scale site-level physiological mechanisms to the landscape scale.

Table 1

Taxon-specific PnET-Succession parameters. Tuned parameters indicated by *. Red maple and red oak were parameterized based on Coulston et al. 2003 because they were not included in the Aspen-FACE experiment

Our study relied heavily on the capabilities of PnET-Succession to simulate the species-specific effects of CO2 and O3 on photosynthetic output. In the PnET-II model, elevated CO2 concentration is assumed to reduce stomatal conductance and increase internal leaf CO2 concentration (Ollinger et al. 2002), thereby reducing transpiration losses and increasing water use efficiency, and photosynthetic capacity (e.g., De Kauwe et al. 2013). This is implemented by computing a CO2 enhancement factor (DelAmax) based on an empirical relationship between the rate of photosynthesis and CO2 concentration. We used a similar approach in PnET-Succession v3.1 to calculate the CO2 enhancement factor, but we used the equation developed by Franks et al. (2013). We slightly modified the Franks’ equation to use internal leaf CO2 concentration (Ci) instead of ambient CO2 concentration (Ca) because Ci better reflects conditions experienced by photosynthetic tissues:

where γ is the CO2 compensation point, assumed to be 40 ppm based on an optimal temperature for photosynthesis of 25 °C, and Cref is 350 ppm.

We replaced the ozone functions of the PnET-II model (Ollinger et al. 1997) to simulate two distinct effects of ozone on photosynthesis: (1) impairment of stomatal function and (2) damage of photosynthetic capacity at the chloroplast level. Elevated ozone induces stomatal sluggishness, preventing complete stomatal closure (Paoletti and Grulke 2005; Hoshika et al. 2012b). Stomatal sluggishness enhances stomatal ozone flux, which causes further damages to photosynthetic tissues. In a previous modeling study, Hoshika et al. (2018b) reported the enhancement of stomatal ozone flux by sluggishness, which could explain the significant ozone-induced decline of photosynthesis in late summer for Japanese white birch. Another factor influenced by stomatal sluggishness is soil water. Hoshika et al. (2015) and Lombardozzi et al. (2015) suggested that sluggishness influences transpiration and thus may change the plant-soil water balance. In fact, several studies (McLaughlin et al. 2007; Sun et al. 2012) reported substantial increases of transpiration at the ecosystem-level in southeastern US forests that were attributed to ozone-induced stomatal sluggishness. This exacerbates the effect of drought because impaired stomata reduce the ability of trees to conserve water (Hoshika et al. 2018a). In the Aspen FACE experiment, the stomatal CO2 response of aspen-birch communities was inhibited by ozone and this effect developed progressively over the growing season (Onandia et al. 2011). As a result, ozone decreases intrinsic water use efficiency and thus increases Ci (e.g., Farage and Long 1995). Using empirical measures of stomatal response to ozone and water stress in a different free-air ozone exposure experiment (Hoshika et al. 2018a; see Paoletti et al. 2017 for the detail of the experiment), we developed equations to modify Ci (Cimodifier) as a function of cumulative ozone dose, water stress and species susceptibility to stomatal sluggishness in response to ozone, to reflect the effects of ozone on stomatal conductance:

where fWater is the reduction multiplier (0.0–1.0) reflecting the water stress being experienced by the cohort, intercept and slope are species-specific intercept and slope coefficients (mol/nmol; Table 3) reflecting each species’ relative susceptibility to ozone-induced stomatal sluggishness derived from stomatal conductance data from an experimental manipulation of ozone in three common European oaks (Quercus ilex, Q. robur and Q. pubescens: Hoshika et al. 2018a), and CumD40O3 is the current month’s cumulative O3 dose above a threshold concentration of 40 nmol/mol, compiled from hourly values (0800–1900 h) since the beginning of the growing season (Ollinger et al. 1997). This modification ripples through the model, reducing the CO2 fertilization effect (DelAmax), reducing water use efficiency and exacerbating water stress, and ultimately, affecting photosynthetic output.

Table 3

Regression coefficients for computing the effects of ozone on stomatal sluggishness (Cimodifier) based on data from Hoshika et al. (2018a, b)

Species-specific stomatal sensitivity to O3

Species used to generate coefficients

Intercept (mol/nmol)

Slope (mol/nmol)

Tolerant

Quercus ilex

0.0087

− 0.021

Intermediate

Q. pubescens

0.0062

− 0.0148

Sensitive

Q. robur

0.0118

− 0.0176

These modifications resulted in explicit tracking of stomatal conductance not found in PnET-II. PnET-Succession now mechanistically simulates the flux of CO2 (JCO2) in and out of the leaves independently from the flux of water (JH2O):

where NetPsn is net photosynthesis (μmolCO2/m2 foliage/s) and Tave (°C) is average monthly temperature (Tmax+Tmin)/2.

Ozone may impair leaf-level physiological functions such as the activity of ribulose-1,5-bisphosphate carboxylase/oxygenase (Rubisco), reducing the photosynthetic capacity (Amax) of the canopy. Evidence from the Aspen-FACE experiment and elsewhere indicates that taxa can vary in their susceptibility to this damage, partly due to antioxidant production (Karnosky et al. 2005). An ozone reduction multiplier (fO3) is applied to reflect this ozone damage by reducing Amax as a function of canopy position (O3 concentration declines down through the canopy profile), length of time of O3 exposure, and stomatal conductance of O3 into the leaves:

where lastO3Effect is O3Effect from the prior month of the current growing season, O3coeff is a species ozone sensitivity coefficient, gwv is stomatal conductance to water vapor, o3 is the growing season cumulative D40O3 including the current month, and relO3 (0.0–1.0) represents O3 extinction through the canopy to the current subcanopy layer. The constant 0.0026 is the kO3Eff coefficient used in PnET-II. The ozone sensitivity coefficient (O3coeff) = 0.0 when there is no ozone-induced tissue damage, and it increases continuously as susceptibility increases, with a value of 3.0 representing three times as much damage as a value of 1.0. PnET-II used a Drought modifier to reflect stomatal closure during drought, but we explicitly simulate stomatal activity (and drought effects) using Cimodifier as described above.

Because of the massive difference in growth rate of clone 216 between dominant and sub-dominant canopy positions in the experiment, we found it necessary to make the primary growth rate parameters dynamic as a function of light availability (Niinemets 1997). First, foliar nitrogen (FolN) determines maximum photosynthetic capacity, and it varies with light reaching each canopy sublayers using:

where AdjFolN is the foliar nitrogen value for the sublayer, MinFolN is the minimum foliar nitrogen value observed for the species, MaxFolN is the maximum foliar nitrogen value expected at full light saturation, fRad is the light stress reduction multiplier (0.0–1.0), and FolNShape is a shape parameter controlling the shape of the FolN adjustment function, allowing both linear and non-linear curves. This function is species-specific to account for differences in physiological and foliar nitrogen response to relative irradiance (Niinemets 1997). Second, foliage fraction (FracFol) determines foliar mass as a fraction of active wood at an annual time step. Similarly, this was also made dynamic as a function of available light such that

where AdjFracFol is the foliage fraction for the current year, MinFracFol is the minimum foliage fraction (under acute light stress), MaxFracFol is the maximum foliage fraction expected at full light saturation, fRad is the light stress reduction factor, and FracFolShape is a shape parameter controlling the shape of the FracFol adjustment function. Finally, we made shade tolerance (HalfSat) dynamic as a linear function of CO2 concentration (Kubiske and Pregitzer 1996) according to:

where AdjHalfSat is the modified shade tolerance parameter, CO2HalfSatEff is a slope coefficient, CO2 is the atmospheric CO2 concentration, and HalfSat is the nominal half saturation light level (units as specified in the climate input file).

Model calibration

We calibrated the model using the results of the Aspen-FACE experiment. Three species (hereafter, taxa) combinations were planted within three subsets (sections) of each treatment ring: (section I) five clones of aspen, one of which (clone 259) died out early in the experiment; (section II) maple and aspen clone-216; and (section III) birch and clone-216. The maple and birch seedlings came from open-pollinated seed sources in northern Wisconsin. There were four treatments, replicated three times: Control, elevated CO2 (+CO2), elevated ozone (+O3), and elevated CO2 and ozone together (+CO2 + O3). CO2 fumigation occurred only during daylight hours, and O3 fumigation followed a prototypical diurnal profile. Experiment details are in Kubiske et al. (2015). Results showed that elevated CO2 tended to enhance growth rates while elevated O3 tended to inhibit growth, reduce vigor and resistance to insects and disease, and increase mortality (Karnosky et al. 2003). However, there were important differences in the response to CO2 and O3 levels among taxa. Some taxa were more O3-tolerant than others, and the response to enrichment was not uniform across taxa (Karnosky et al. 2005). Furthermore, various clones of aspen had divergent responses to the treatments, including one (clone 8) that fared better under elevated O3 and worse under elevated CO2 because its competitors had a greater response to those treatments (Kubiske et al. 2007; Moran and Kubiske 2013).

We calibrated PnET-Succession using the individual treatment factors and validated it against the combined factor treatment (+CO2+O3). Treatments were simulated by modifying CO2 and O3 inputs in the monthly weather time series and keeping all other parameters and inputs unchanged. We simulated each Aspen-FACE replicate using the CO2 and O3 concentrations measured within its experimental plot. We also assigned each replicate number (N = 3) a specific random number seed to allow comparison of model variability relative to replicate treatment variability given that there was no treatment variability among the Control replicates. Hourly CO2 values ranged between 350 and 390 ppm on ambient CO2 plots and 502–530 ppm on elevated CO2 plots (Kubiske et al. 2015). Hourly O3 measurements ranged from 32 to 38 ppb on ambient O3 plots to 38–49 ppb on elevated O3 plots.

Initial calibration of PnET-Succession was conducted by comparing model output of each control Aspen-FACE replicate treatment to the empirical measurements made in the corresponding replicate control plot. We simulated the years of the experiment (1998–2008) plus 5 years, with a site (cell) size of 30 m, allowing simulation of an experimental plot with a single cell. Each cell was initialized with 1-year-old cohorts of the taxa combinations of the Aspen-FACE experiment. Clone 259 was not included because there were insufficient data to estimate its model parameters. We used weather data from a National Weather Service observation station located within 16 km of the experimental site (Rhinelander WI airport). Model parameters were initially set to their empirically derived values as published for other studies (e.g., Aber et al. 1995, 1996; Gustafson et al. 2017), and then calibrated to match Aspen-FACE observations of aboveground live biomass, foliage biomass, and leaf area index (LAI) within the 95% confidence intervals of empirical measures from each Aspen-FACE experimental plot (Kubiske et al. 2007; Kubiske 2013). We iteratively modified (tuned) eight species parameters (Table 1) to produce behavior consistent (both in magnitude and temporal variation) with the empirical measurements under ambient (control) conditions. Because foliar nitrogen (FolN) is the primary determinant of growth rates in PnET-Succession (Aber et al. 1995, McKenzie et al.: Local and global parameter sensitivity within an ecophysiologically based forest landscape model, in review), we began with minimum and maximum foliar nitrogen values found in the literature for each species and tuned FracFol values to produce the observed amounts of foliage for each taxa (for acronym definitions see Table 1). We then tuned FolN and FracActWd to produce the observed amount for woody and foliar biomass for each taxa to reflect the observed differences in growth rates among taxa, and then modified SLWmax (within empirical limits found in the literature) to mimic the observed LAI as closely as possible, given the foliage mass present. The site parameter MaxDevLyrAv controls the assignment of cohorts to canopy layers based on differences in their size (biomass), and it was calibrated to achieve canopy layering in the year in which it was observed in the Aspen-FACE experiment. We strove to use the same or similar values for each parameter across all species where possible, to mitigate the confounding effects of parameter variation on the competitive interactions between taxa.

We then used the model to simulate the individual treatment factors. CO2 response is primarily controlled by the modified equation (Eq. 1) of Franks et al. (2013). The CO2HalfSatEff parameter was set a priori by shade tolerance (HalfSat) class, with less effect of elevated CO2 on HalfSat for more shade tolerant taxa. The O3 response was calibrated by tuning the O3StomataSens and O3GrowthSens input parameters (Table 1) to mimic the observations on the +O3 treatment replicate plots.

Model validation

After completing calibration, we validated model behavior by simulating the +CO2+O3 treatment and comparing model output to empirical measurements of that treatment. Similar to model calibration, empirical measures from the treatment plots were compared to those output by the model, and acceptable agreement was assumed when the predicted values passed within or near the 95% confidence intervals of the time series of empirical measurements, recognizing that the model was calibrated to increase generality for landscape modeling purposes rather than to predict the results of this single experiment. We also visually assessed whether the relative magnitude and trend of a predicted variable was consistent with observed patterns.

Extrapolations of the Aspen-FACE experiment

Using the validated model, we virtually extended the Aspen-FACE experiment 80 years into the future for the forest communities used in the experiment and expanded the scope of the experiment by simulating all six taxa together. We also generated a hypothetical assemblage of species that might co-occur with those used in the Aspen-FACE experiment (Table 1), and grew them together for 150 years. Ozone sensitivities of species not included in the Aspen-FACE experiment were parameterized based on classifications in Coulston et al. (2003). In each virtual experiment, all taxa had an initial age of 1 year and were exposed to the same site and weather and atmospheric conditions used for calibration and testing. The weather and atmosphere time series for each treatment combination were repeated to create a 150-year time series, replicating the experiment weather for the entire 150 years (i.e., no climate change). Each experiment was conducted on a single-cell with dispersal, establishment, and disturbance turned off to eliminate their confounding effects. Such effects will be investigated in a forthcoming study, and our results should be interpreted in light of the fact that there was no establishment of new cohorts.

Analysis procedures

Quantitative testing of model performance was problematic because the model predicts growth trends at a monthly time step and the empirical measures were estimates based on samples of physical traits taken at various times during a growing season. We therefore relied on visual evaluation of graphical overlays of simulated trends and empirical measures, rather than conducting statistical tests of our hypotheses, as recommended for simulation experiments by White et al. (2014), concentrating on the relative magnitude of predicted values and measurements and their variation through time. This approach allowed us (and readers) to visually assess how consistent the simulation results were with observations. Variables evaluated were above ground woody biomass, foliage biomass, and LAI, and we plotted the observed value for each replicate and the 95% confidence interval for comparison with simulated values.

Model calibration

Calibrating individual taxa under the control treatment was fairly straightforward with the tuning parameters used, allowing us to almost always produce model predictions that fell within the 95% confidence intervals of empirical measures of taxon growth (Fig. 1). Foliage biomass sometimes differed from empirical observations because foliage and growth rates were not always consistently related in the empirical observations. LAI values were computed from foliage biomass values. Clone 216 was interesting because it occurred in all three taxon combinations, and its behavior was quite different in each, depending on its competitors. We were therefore forced to calibrate clone 216 such that its parameters represented a compromise of its behavior with different co-occurring taxa.

The +CO2 treatment increased each response variable by an amount comparable to that observed in the empirical experiment, producing an agreement with empirical measures very similar to that seen under the control treatment (Fig. 2). Individual taxa were not always assigned to the same canopy layer in each replicate, which sometimes produced a divergence among replicates. We were able to calibrate the model to mimic the observed decline of the aspens (section I) to the +O3 treatment, including the release of the ozone tolerant clone 8 (Fig. 3). The relative response to ozone in the other sections was appropriately predicted, but the absolute response was less satisfying.

Model validation

The +CO2+O3 treatment constituted an independent test (validation) of the calibrated model (i.e., there was no calibration to the +CO2+O3 treatment), and the predictions generally fell within the 95% confidence intervals of observations and the relative differences in woody biomass of the simulated taxa conformed to the observations (Fig. 4). Foliage biomass and LAI followed woody biomass as they did for the other treatments. Again, the relative responses to the treatment were appropriately predicted, but the absolute responses were sometimes less satisfying, not always passing through the 95% confidence intervals.

Extrapolations of the Aspen-FACE experiment

Extrapolation of the Aspen-FACE experiment 80 years into the future showed that in the aspen section, clone 271 continued its aggressive growth over the longer time period, and that the other taxa thrived or languished depending on treatment, with clone 8 being slightly more dominant than clone 216 under treatments with elevated ozone (Fig. 5). In the birch and maple sections, the treatments affected the relative ability of taxa to grow, but did not change the relative dominance of any taxon (Figs. 6 and 7), although in a single replicate, clone 216 was dominant over birch under the +CO2+O3 treatment.

Fig. 5

Extrapolation by simulation 80 years into the future of the aspen section of the Aspen-FACE experiment. Error bars show one standard deviation of three replicates

Fig. 6

Extrapolation by simulation 80 years into the future of the birch-aspen Section of the Aspen-FACE experiment. Error bars show one standard deviation of three replicates

Fig. 7

Extrapolation by simulation 80 years into the future of the maple-aspen Section of the Aspen-FACE experiment. Error bars show one standard deviation of three replicates

When we put all six taxa together, the competitive outcomes varied dramatically by treatment (Fig. 8). Clone 271 outcompeted birch under the control treatment, but three clones outcompeted birch under elevated CO2. Birch and clone 8 dominated under the +O3 treatment as expected given their relative tolerance of O3. Maple was slow-growing and very subdominant in the Aspen-FACE experiment, and our extrapolation indicates that its growth was just adequate to survive under all treatments, although it survived in only one replicate under the +O3 treatment.

Fig. 8

Extrapolation by simulation 80 years into the future of all six Aspen-FACE experiment taxa together. Error bars show one standard deviation of three replicates, and a taxon survived in only one replicate where there are no error bars

The hypothetical assemblage caused a novel collection of species to grow and compete under the experimental treatments, and the results demonstrated that outcomes varied widely depending both on treatments and stochasticity of competitive interactions (Fig. 9). Where a mean trend line falls within the error bars of another taxon, this indicates that these taxa sometimes traded dominance in those time steps. Note that mortality occurring within approximately 30 years of longevity age is at least partially the result of senescence, while mortality occurring prior to that is the result of competition failure. Generally, oak was dominant, but the treatments modified the ability of each taxon to survive and compete, with elevated CO2 boosting the productivity of fast-growing pioneer taxa (e.g., clone 271, balsam poplar).

Fig. 9

Extrapolation by simulation 150 years into the future of an arbitrary assemblage of species found in northern Wisconsin under the treatments of the Aspen-FACE experiment. Error bars show one standard deviation of three replicates, and a taxon survived in only one replicate where there are no error bars

We identified a number of areas in which our model achieved our objectives. (1) We were able to tune a small number of species parameters to produce woody biomass within or near empirically observed 95% confidence intervals for all taxa under the control treatment. (2) The Franks et al. (2013) CO2 effect equation in the model (coupled with the modest CO2 effect on shade tolerance) was able to mimic the empirical CO2enhancement for the CO2 concentrations seen in the experiment. (3) Similarly, the response of all taxa to the +O3 treatment mimicked empirically observed declines in both relative and absolute terms. (4) The model successfully predicted the outcome of the +CO2+O3 treatment, which served as an independent validation of the model. (5) Predictions of relative treatment effects were very consistent with those seen in the Aspen-FACE experiment (Table 3). (6) The model extrapolated treatment effects for a number of taxon combinations and provided insight into the relative long-term effects of the treatments and stochasticity on competitive outcomes. (7) We were able to simulate CO2 and O3 effects without a large computational overhead, retaining consistency with the model’s purpose for landscape-scale simulations. This suggests that our modifications to PnET-Succession now enable it to account for a broader range of stressors (O3 pollution) in projections of future forest dynamics for research and management questions.

We also identified areas where our model failed to meet our objectives or did not accurately predict all cohort state variables, and these should be understood to properly interpret our results and prior to applying the model to predict landscape-scale responses to elevated CO2 and O3. (1) Although we modified the model to make growth parameters dynamic in response to light environment, we were still unable to find a single set of input parameters for clone 216 that precisely mimicked empirical growth observations in the presence of all competitors under all treatment combinations. The discrepancy was most pronounced for clone 216 growing with maple under the +CO2 treatment, although it should be noted that the empirical measures were incredibly high for a 10-year-old cohort. (2) The model did not always reliably predict foliar biomass (and therefore, LAI). Here we note that the relationship between foliage biomass and biomass growth was not consistent among taxa and treatments, making it difficult to model. PnET-Succession allocates a deterministic portion of net photosynthates to foliage, computed as a fraction of the active wood biomass that delivers water to the canopy. This apparently produced a reasonable estimate of foliar mass for most of the taxa combinations, but the agreement was inconsistent. (3) When model variability was controlled (constant random number seed) and each simulation replicate used the concentrations of CO2 and O3 observed in the experimental replicates, replicate variability was only about twice the width of the trend lines in Figs. 1, 2, 3, and 4 (not shown). However, our model did not include any of the uncontrolled factors likely to exist in the empirical experiment such as belowground abiotic and biotic conditions, plot scale disturbances (e.g., pests, herbivores), and genetic variation of planted stock. For example, ozone may alter the feeding behavior of insects, affecting foliar biomass (Agathokleous et al. 2017). Our results suggest that these factors account for most of the empirical replicate variability rather than variability in the application of the treatments.

Prior versions of PnET-Succession assumed that foliage biomass is proportional to active woody biomass (transport tissue) rather than light availability. We discovered that this was inadequate to simulate the FACE experiment because growth patterns varied greatly among aspen clones and between clone 216 and paper birch (Moran and Kubiske 2013). For example, clone 216 exhibited increasing incidence of determinate growth through time, and almost no sylleptic shoots (shoots that form opportunistically throughout the growing season) formed later in the experiment. Likewise, aspen clone 259, which died out after 6 years, was greatly suppressed and exhibited determinate primary growth (M.E. Kubiske, personal observations). In contrast, clone 271 produced copious indeterminate, sylleptic shoots throughout the experiment, and eventually became the dominant clone in all treatments. Similarly, paper birch retained indeterminate primary growth throughout the experiment and also tended to produce sylleptic shoots. We speculate that these growth patterns of clone 271 and birch were the direct cause of their much greater foliar biomass and rapid accumulation of biomass. Conversely, juvenile sugar maple, with determinate primary growth, is well known to have very slow height growth for several years, after which its determinate growth rate increases (Burns and Honkala 1990); thus, clone 216 quickly over-topped sugar maple early in the experiment and maintained dominance. The original PnET-II model dynamically allocated foliage to successively lower canopy strata until light availability is inadequate to support photosynthesis. This approach was not included in PnET-Succession because of the computational overhead of this approach for the millions of cohorts (each with multiple canopy sub-layers) that often comprise a landscape simulation. Although our modifications to the model did not allow us to perfectly mimic empirical observations of the Aspen-FACE experiment, they did produce credible behavior using simplified functions that allow the model to be tractable at landscape scale.

Our results demonstrate that long-term competitive outcomes can vary dramatically depending on atmospheric conditions, consistent with a 2-year phytotron study showing that ozone sensitivity of trees may affect competitiveness of species depending on their ability to enlarge crown volume (Kozovits et al. 2005). Perhaps, the most dramatic example in our 150-year simulations was the performance of red oak (Fig. 9). Red oak was a dominant presence in all the treatment combinations except under elevated CO2, where it was subordinate to the fast-growing clone 271 (Fig. 9). The tolerance of clone 8 to elevated ozone allowed it to remain completive for 80 years in all scenarios that included ozone, except when competing with red oak, which is also ozone tolerant, and also more shade and drought tolerant. Similarly, in the simulation of all Aspen-FACE taxa together, the somewhat ozone tolerant birch fared best under treatments with ozone, usually at the expense of clone 42 (Fig. 8). The model also makes some non-intuitive predictions. Although birch was dominant over clone 216 in the 11-year Aspen-FACE experiment and when extrapolated for 80 years, when competing with other taxa, birch was not able to outcompete the less shade tolerant clone 271 except in the presence of elevated ozone (Fig. 8). Birch was also completely unable to compete with red oak (Fig. 9), likely because of oak’s greater shade, drought, and ozone tolerance.

Caveats

Our model simulates elevated tropospheric ozone as a chronic disturbance at a monthly resolution, as opposed to occasional spikes in O3 concentration that commonly reach forested areas with a daily or weekly resolution. The +O3 treatment of the Aspen-FACE experiment had a target concentration of 1.5× ambient (generally 45–65 ppb O3), and occasional spikes of much higher O3 naturally impacted the site. In addition, no O3 was dispensed into the +O3 treatment plots on heavily overcast or rainy days, when air temperature was below 15 °C, or when wind speeds exceeded 0.4 m s−1. Thus, the actual concentrations in the +O3 treatment, integrated over the course of a growing season, were within 10% of target (above or below) 66% of the time and within 20% of target 83% of the time. PnET-Succession is unable to simulate O3 spikes at the weekly time scale nor does it vary O3 with weather conditions. We hypothesized that delayed canopy closure due to elevated O3 could result in recruitment of greater understory biomass as observed in the Aspen-FACE experiment, although we did not simulate this process.

The mechanisms of the effects of elevated CO2 on photosynthesis in PnET-Succession are fairly simplistic. In the model, elevated CO2 increases internal CO2 concentration (Ci) and reduces stomatal water loss, increasing water use efficiency with increasing CO2. Elevated CO2 also modestly increases shade tolerance (Kubiske and Pregitzer 1996), which has an interacting effect on FolN and FracFol, which are dynamic as a function of light stress in PnET-Succession v3.1. We experimented with a direct link between CO2 and FolN, but were unable to find a scientifically defensible algorithm to implement such an effect short of adding full N cycling to the model.

The Aspen-FACE experiment was just one replicate at landscape scale, so the parameterization of our model based on those results should be tested elsewhere. However, because data sets for such tests are exceedingly rare (Karnosky et al. 2007), our results provide an important launching pad for landscape simulation of the effects of atmospheric pollution on forests. PnET-Succession now provides a powerful tool for scaling such plot-level empirical results to landscape scale. While our estimates of O3 response parameters for taxa included in the experiment are fairly certain, those parameters for species not included in the experiment are less certain, being based solely on generic estimates of ozone tolerance published by Coulston et al. (2003). Nevertheless, these estimates enable a mechanistic and proportional response to ozone that allows study of the generic response of forest communities to ozone pollution at landscape spatial and temporal scales, a capability that has hitherto been lacking.

Our results demonstrate the power of a mechanistically motivated approach to predict the competitive outcomes of species with varying life history characteristics under novel atmospheric conditions for long time periods. Growth rate is a key driver of the competitive interactions among taxa, being the primary determinant of the dynamics of young forests. However, over longer time frames (centuries), species longevity and shade tolerance can change forest composition in profound ways that also interact with disturbance regimes (Gustafson et al. 2017). Also, taxa do not always perform over the long term as they do over the short term (e.g., compare clone 216 (growing with other aspen clones) under +O3 between Figs. 3 and 5). Furthermore, it is very difficult to make robust long-term predictions of tree species competitive outcomes under novel abiotic conditions such as climate change and elevated atmospheric pollution, or novel assemblages of species such as seen in Figs. 8 and 9. In our study, we were able to use a mechanistically-motivated landscape model to make predictions that were consistent with results from the short-term Aspen-FACE experiment (Kubiske et al. 2007; Moran and Kubiske 2013). Our use of the validated model to extrapolate to longer time frames and to novel assemblages of taxa provides an important demonstration of the use of mechanistic models to scale plot-level experimental results to longer time frames and to other species not included in the empirical experiments. Applications to real landscapes will present new challenges, such as choosing a single set of parameters for a clonal species like aspen that exhibit great variation in some parameters, such as ozone sensitivity.

Simulating ozone impacts on forest competition is challenging, and this is a major contribution of our study. Our study is quite similar in design to a previous study (Gustafson et al. 2013) that used the more phenomenological LANDIS-II succession extension, Biomass Succession (Scheller et al. 2007), and it is instructive to compare these studies. Biomass Succession constrains maximum biomass for each taxon based on generic estimates for average abiotic conditions and does not simulate competition for light and water for growth, instead modeling a generic competition for growing space. The present study predicted much higher long-term biomass of taxa than did our previous study because we used a model where biomass is an emergent property of physiological attributes, monthly abiotic inputs, and mechanistic simulation of photosynthesis. Our prior study did not include an independent test of the model because all four treatment combinations had to be calibrated. In both studies the growth potential of taxa was the primary determinant of short-term taxon dominance, but longevity, shade, and drought tolerance had an important effect in the longer term. In our prior study, the long-term difference between clone 216 and birch or maple was extremely high, reflecting the large differences seen in the short-term because Biomass Succession does not account for the competition for light and water for growth, instead modeling only competition for growing space. In both studies, taxa responded differently to each treatment. However, no taxa were in our prior study except by senescence, even under treatments for which they were poorly adapted, unlike in the present study. When all 6 taxa were grown together in our prior study, the general trends were similar to the present study for the +O3 and +CO2+O3 treatments, but for the control and +CO2 treatments, clone 217 was even more dominant, again because of differences in how competition was simulated. These differences underscore the value of using a more mechanistic approach that better simulates the processes that determine competitive outcomes at the appropriate temporal resolution of the drivers.

Our modifications to the model included the effects of CO2 and O3 on photosynthesis by explicitly modeling the effect of O3 on Ci, and the model produced good agreement with observed biomass trends in the Aspen FACE experiment. Most previous modeling studies assumed a tight coupling of stomatal conductance and photosynthesis (e.g., Martin et al. 2001), because stomatal conductance is generally regulated so as to maintain a fixed ratio of Ci to ambient CO2 concentration (Larcher 2003). However, an increase of Ci was often observed after O3 exposure in experimental studies (Hoshika et al. 2013, 2014; Watanabe et al. 2014), and this has been attributed to ozone-induced stomatal sluggishness (Paoletti 2005; Mills et al. 2009; Hoshika et al. 2014). This ozone-induced sluggish response to CO2 was also observed in the Aspen FACE experiment (Onandia et al. 2011). A recent study found that the expression of genes that signal the stomata to close in response to elevated CO2 is inhibited by O3 (Dumont et al. 2014). Watanabe et al. (2014) reported a decrease in the ratio of photosynthesis to stomatal conductance for Monarch birch (Betula maximowicziana) exposed to O3, in tandem with an increase in Ci resulting from relatively high stomatal conductance, which might be a compensation response against an O3-induced decline of photosynthetic capacity in chloroplasts. Although the underlying mechanisms of stomatal sluggishness are still under investigation (e.g., Wilkinson and Davies 2010), the addition of stomatal sluggishness effects on Ci into models decouples photosynthesis and stomatal conductance under elevated O3 and CO2 conditions, with potentially large impacts on carbon and water fluxes at the landscape level.

Our primary purpose in this study was to improve and test the ability of a mechanistically motivated forest landscape model to predict the growth and competitive interactions of tree species to the globally changing drivers of CO2 concentration and ozone pollution, with the ultimate goal being the study of how these pollutants might interact with disturbances and other spatial processes that structure forests to impact forest composition and spatial pattern at landscape scales. We therefore strove to not over-parameterize the model to precisely fit a single set of empirical data, but to use parameter values that were similar among taxa so that they could readily be applied on heterogeneous landscapes under a wide range of atmospheric conditions. Additionally, we did not wish to encumber the model with so much computational burden that it would be infeasible to run at landscape scales over prediction horizons of centuries. One consequence of this approach is that the ability of the model to match empirical observations for all state variables was not as strong as it might otherwise be. Nevertheless, the prediction of the cohort attribute most important for landscape simulations of growth and competition (woody biomass) was generally quite good, and the relative responses (among taxa) were quite consistent with empirical observations even when the absolute responses were not exactly consistent. In all cases, the competitive outcome of all taxa combinations under all treatments was identical to that seen in the empirical observations. We conclude that PnET-Succession v3.0 can generate potentially robust predictions of the effects of elevated CO2 and O3 at landscape scales because of its mechanistic algorithms, but this should be confirmed by tests against other data sets.

Acknowledgements

We thank Scott Ollinger and Zaixing Zhou for helping us understand the ozone functions in PnET-II. Thanks to Brian Sturtevant, Andy Burton, Scott Ollinger, and two anonymous reviewers for critique of the manuscript.

Funding

Funding was provided by the Northern Research Station of the USDA Forest Service. The Aspen-FACE experiment was principally supported by the Office of Science (BER), US Department of Energy Grant No. DE-FG02-95ER62125 to Michigan Technological University; Contract No. DE-AC02-98CH10886 to Brookhaven National Laboratory; Office of Science (BER), US Department of Energy Interagency Agreement No. DE-AI02-09ER64717 to the US Forest Service, Northern Research Station; the US Forest Service Northern Global Change Program; and the Canadian Forest Service.

Availability of data and materials

The datasets supporting the conclusions of this article are included within the article (and its additional files). The model is freely available at http://www.landis-ii.org.

Authors’ contributions

EJG conceived the study and drafted the manuscript. MEK designed the algorithms for model modification. BRM programmed the model. YH and EP provided data and made writing contributions. All authors read and approved the final manuscript.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Gustafson EJ (2013) When relationships estimated in the past cannot be used to predict the future: using mechanistic models to predict landscape ecological dynamics in a changing world. Landsc Ecol 28:1429–1437View ArticleGoogle Scholar

McLaughlin SB, Nosal M, Wullschleger SD, Sun G (2007) Interactive effects of ozone and climate on tree growth and water use in a southern Appalachian forest in the USA. New Phytol 174:109–124View ArticleGoogle Scholar