Collaborative Research: Thresholds and mechanisms of net ecosystem production (NEP) resilience following moderate disturbance: Why does one ecosystem recover and another one crash? (Award Abstract #1655095)

I. Project Description & Background: Aging Forests, Shifting Disturbance Regimes

The structure and function of forests in the US upper Midwest and northeast were dramatically altered by stand-replacing disturbance over a century ago (Frelich and Reich 1995). The resulting cohort of early successional trees covers an immense area (>400,000 km2 of the contiguous US, USDA 2011), and is beginning to decline, while longer-lived species representation is increasing (Hill et al. 2005, Bergen and Dronova 2007, Gough et al. 2010). At the same time, forest disturbances in the region have largely transitioned away from severe events that cause stand replacement, to more moderate disturbances resulting in only partial canopy defoliation or species loss (Pregitzer and Euskirchen 2004, Birdsey et al. 2006b, Luyssaert et al. 2008, Radeloff et al. 2012, Williams et al. 2012, Canham et al. 2013, Lu et al. 2015, Cohen et al. 2016). These subtler disturbances include partial harvests, wind, insects, pathogens, and age-related senescence (Birdsey et al. 2006a, Clark et al. 2010, Hicke et al. 2012, Williams et al. 2012, Cohen et al. 2016), which contribute to a gradient of disturbance intensities across the landscape. Unlike stand-replacing disturbance, long a focus of ecologists (e.g., Kira and Shedei 1967, Odum 1969, Bormann and Likens 1979), less studied moderate disturbances tend to increase canopy structural complexity (Hardiman et al. 2011), with entirely different functional consequences for ecosystems (Peters et al. 2013, Dietze and Matthes 2014, Anderegg et al. 2015, Trumbore et al. 2015).

Although these successional dynamics are well understood, aging and moderately disturbed forests contain surprises for ecologists. In many forests, production shows unexpected resilience to moderate severity disturbance (Amiro et al. 2010, Luo and Weng 2011, Coomes et al. 2012, Hicke et al. 2012, Gough et al. 2013, Peters et al. 2013, Vanderwel et al. 2013, Flower and Gonzalez-Meler 2015). Net primary production (NPP), and consequently net ecosystem production (NEP), may be sustained or increase in older forests that experience moderate disturbance (Gough et al. 2013, Bond-Lamberty et al. 2015), a finding consistent with the intermediate disturbance hypothesis (Huston 2014) and classic successional theory (e.g., Odum 1969) but at odds with predictions of declining production in aging forests (Birdsey et al. 2006b, Wear and Coulston 2015). Long-term time series (Keenan et al. 2013b) from century-old forests show production is increasing, while recent syntheses, including one co-authored by the PIs, found little or no decline in temperate forest production with age or moderate disturbance (Gough et al. 2016, Hudiburg et al. 2009, Amiro et al. 2010, He et al. 2012, Stephenson et al. 2014). Higher-thanexpected production in aging forests is observed in temperate deciduous (Hardiman et al. 2013b, Gough et al. 2016) and coniferous forests (Hudiburg et al. 2009, Meigs et al. 2009, Lutz et al. 2014, Harmon and Pabst 2015), with disturbance severity hypothesized, but not demonstrated, to differentiate C cycling trajectories (Urbanski et al. 2007, Meigs et al. 2009, Gough et al. 2016). These recent observations of sustained or minimally declining NEP in old temperate forests contradict influential age-NEP syntheses (Pregitzer and Euskirchen 2004, Luyssaert et al. 2008) that imply universally large (> 50 %) age-related declines in NEP as a result of lower NPP (rather than > respiration, Tang et al. 2014, Zhou et al. 2015); this highly cited generalized pattern of age-NEP-NPP derives from analyses of temperate deciduous and coniferous forests and (in Luyssaert et al. 2008) biomes, with agerelated production declines driven by >100-yr-old coniferous forests more likely to experience high severity disturbance (Williams et al. 2014, Cohen et al. 2016; and Pregitzer and Euskirchen 2004, Luyssaert et al. 2008 supplementary material). The assumption of large universal declines in NPP and NEP with age features prominently in standard ecology textbooks (Chapin et al. 2002), and has been included in policy documents intended to inform ecosystem carbon management and land-use policy (e.g., see Birdsey et al. 2007; State of the Carbon Cycle Report).

Recent results from the Forest Accelerated Succession Experiment (FASET) at the University of Michigan Biological Station (UMBS) also point to unexpected resilience of aging, moderately disturbed forests (Cote and Darling 2010). This experiment, in which all early successional canopy trees were killed within a 39 ha area, reported a 3-4 year, 44% decline in leaf area, but no corresponding change in NEP, because NPP remained stable despite multi-year partial defoliation (Gough et al. 2013); in contrast, a series of models predicted large short-term declines in NEP owing to a reduction in NPP rather than ecosystem respiration (Bond-Lamberty et al. 2015). Using FASET’s spatial (plot-to-plot) variation in tree mortality Stuart-Haentjens et al. (2015) found that wood NPP exhibited surprisingly high resistance to increasing disturbance severity until a threshold of ~60% basal area mortality was exceeded (Fig. 1). Hardiman et al. (2011, 2013b), examining sites at UMBS, proposed a mechanism for such resilience following moderate disturbance severity, showing that rising canopy structural complexity with increasing forest age corresponded with higher NPP, concluding that the emergence of more complex, older forests may, against expectations, increase their C sink strength. Hardiman (2013) and Stuart-Haentjens (2015) reasoned that rising structural complexity may improve ecosystem resourceuse efficiency, thereby offsetting – to a point – a reduction in canopy C fixation, and consequently NPP, associated with reduced LAI.

Figure 1. The response of wood NPP to increasing disturbance severity at UMBS. Plots are nested within the 39 ha FASET site. A breakpoint of ~0.6 (P < 0.05) suggests a change in the responsiveness of wood NPP to rising disturbance severity. Mechanisms underlying resilience or decline are uncertain. Stuart-Haentjens, Gough et al. 2015, Ecology.

The poor model performance documented by Bond-Lamberty et al. (2015) is not unexpected. Many ecosystem- to global-scale models, designed for and tested in early- to mid-successional forests with low biological and structural complexity and low disturbance frequency and severity, have trouble reproducing late-successional forest C cycling dynamics (Landsberg and Waring 1997, Raulier et al. 1999, Paustian et al. 2000, Law et al. 2003, Li et al. 2003, Zhao et al. 2009, Hicke et al. 2012, Dietze and Matthes 2014). For example, models such as PnET (Peters et al. 2013), ED, ZELIG, and Biome-BGC all predict significant short-term drops in production with moderate disturbance (Bond-Lamberty et al. 2015), in stark contrast to real-world observations. For this reason several recent syntheses have called for replicated, manipulative experiments closely linked with modeling analyses to investigate responses of C fluxes to different disturbance severity levels, and identify specific mechanisms supporting forest resilience (Amiro et al. 2010, Goetz et al. 2012, Hicke et al. 2012, Becknell et al. 2015, Flower and Gonzalez-Meler 2015, Frank et al. 2015, Medlyn et al. 2015, Reyer et al. 2015).

These results highlight significant knowledge and data gaps: What mechanisms determine the threshold between NPP resilience and decline? What are the implications for models? Do temperate deciduous and coniferous forests, with different disturbance regimes, follow unique age-production (NPP and NEP) trajectories? And what data are most needed to resolve these questions? We propose a tripartite approach. Coupled field and modeling experiments will manipulate disturbance severity to identify the mechanisms responsible for this apparent resilience, and their limits, while using an NSF-funded (Award #1458021) data assimilation framework to test how well different models reproduce observed behavior, and why. This goes well beyond previous observational studies or, as in FASET, experiments that applied only a single level of moderate disturbance, often without replication or controls (Amiro et al 2010). Finally, we will perform an ‘open data’, fully reproducible synthesis of North American temperate forests to lay the foundation for future work, by our team and others, examining how and why successional trajectories of NPP and NEP–posited by ecological theory (e.g., Odum 1969) and observations (e.g., Williams et al. 2014, Gough et al. 2016) to be closely coupled with disturbance frequency and severity–differ among forest types. This synthesis will focus on middle- to late-successional forests since this region of the age-production curve is the least studied and understood.

Figure 2. Rugosity, a metric of canopy structural complexity, is positively correlated with ANPP (a) and improved light-use efficiency (LUE, b). A similar relationship (not shown) was found between rugosity and nitrogen-use efficiency. More structurally complex and resource-use efficient forests that emerge from moderate disturbance may thus resist declines in production up to some critical threshold of disturbance. Hardiman et al. 2013, For. Ecol. Man.

II. Objectives and Hypotheses

O1: Identify mechanisms behind both NPP resilience and decline following disturbance.

  • H1.1: Resilience of NPP to disturbance will be caused by canopy structural changes that increase resource-use efficiency; declines will occur when improved resource-use efficiency is outweighed by severe reductions in LAI and thus canopy photosynthesis.

    Gough et al. (2013) observed a disturbance-related reduction in landscape-scale LAI of nearly half with no change in NPP. This high resistance of NPP to LAI loss corresponded with improved light-use efficiency and sustained high rates of canopy light absorption. Observational work conducted at UMBS across undisturbed plots varying in canopy structural complexity indicates that more multilayered, heterogeneous (i.e., rugose) canopies have higher light and nitrogen use efficiencies, and are thus more productive (Fig. 2, Hardiman et al. 2013). Such structurally complex canopies may use light more efficiently than simpler canopies, contributing to higher NPP (Ahl et al. 2004, Martin and Jokela 2004, Walcroft et al. 2005, Duursma and Makela 2007, Leuschner et al. 2009, Niinemets 2010, Stark et al. 2012). However, the quantity and surface area of foliage is an essential predictor of production in many ecosystems (van Dijk et al. 2005, Luyssaert et al. 2007, Lindroth et al. 2008, Duursma et al. 2009, Reich 2012), suggesting that critical non-linear thresholds exist (Goetz et al. 2012) in which partial canopy defoliation cannot be offset by improved resource-use efficiency.

  • H1.2: Improvements in whole-canopy light-use efficiency following disturbance will be driven by rapid leaf-level physiological and morphological changes that permit canopy subdominants to take advantage of new light.

    Canopy structure affects light transmission through, and interception by, the canopy (Parker and Russ 2004); consequently, forests with greater structural complexity (i.e., multi-layered with more gaps) are less susceptible to self-shading, improving light harvesting by subdominant canopy elements and thus light-use efficiency of the whole canopy (Gough et al. 2007a, Campbell et al. 2009). Such forests distribute light more evenly throughout the canopy (Walcroft et al. 2005, Duursma and Makela 2007, Pangle et al. 2009), improving mean leaf-level photosynthesis because of the non-linear relationship between C assimilation and absorbed photosynthetically active radiation (Parker et al. 2001, Parker and Russ 2004, Walcroft et al. 2005, Bartemucci et al. 2006, Niinemets 2007, Stark et al. 2012, Niinemets et al. 2015). As canopy gaps form, less solar radiation is available to drive upper canopy photosynthesis, but with minimal resulting declines in mean leaf-level C assimilation by this light-saturated stratum (Canham et al. 1999, Ishii and Asano 2010). Improved light availability and physiological functioning of the lower canopy thus may offset declines in the photosynthetic contribution of upper canopy leaves following a moderate reduction in leaf area, but this is a particularly difficult dynamic for most fixed-parameter models to reproduce (Dietze 2014). Lower canopy plants are important contributors to C storage in forests, including those recently subjected to upper canopy defoliation (Ueyama et al. 2006, Gough et al. 2007a, Misson et al. 2007, Campbell et al. 2009, Stuart-Haentjens et al. 2015). This is a putative, and presently unconfirmed, mechanism for sustained and perhaps higher rates of NPP in old forests (Hardiman et al. 2011, 2013).

  • H1.3: NPP resilience to disturbance will be greatest when nitrogen (N) from senescent individuals is reallocated in support of LAI recovery, and moderately disturbed canopies emerge more structurally complex than undisturbed forests.

    LAI recovery is an important process supporting production resistance and resilience to disturbance in forests (Amiro et al. 2010), with the replacement of leaf area, and thus photosynthetic capacity, lost to disturbance closely coupled with canopy N reallocation (Nave et al. 2011). In the first phase of the FASET experiment, leaf area in the treatment forest recovered to predisturbance levels as N from senescent aspen and birch was fully reallocated, rather than leached or immobilized, to residual canopy species in support of new foliage production (Nave et al. 2011, Gough et al. 2013, Nave et al. 2014). Studies at UMBS show that NPP is positively correlated with both canopy N mass and structural complexity, with greater structural complexity and canopy N mass together conferring high rates of production (Hardiman et al. 2013b). Importantly, several recent studies suggest physical and biological complexity confer functional resilience in forest ecosystems, but without explicitly identifying underlying mechanisms (Seidl et al. 2014, Pedro et al. 2015, Metz et al. 2016).

O2: Understand why models fail to replicate observed NPP resilience.

  • H2.1: Neither a classic big-leaf model (Biome-BGC) nor a partial-differential gap model (ED2, as well as the R version of ED incorporating ecophysiological disturbances; Dietze and Matthes, 2014) will replicate observed NPP resilience across a disturbance severity gradient.

    Modeling studies generally either parameterize (Weng et al. 2012) or produce (Peters et al. 2013) mildly nonlinear responses, but exhibit little resilience to moderate disturbance (Dietze and Matthes 2014). We expect that Biome-BGC will be improved by incorporating metrics of canopy structure (proxies for aggregated changes in leaf physiology and shifts in whole canopy light-use efficiency, H2.1), while ED’s lack of a nitrogen cycle will compromise its performance. Following H2.1 and 2.2 above and our previous work in this area (Fig. 3), we anticipate that the treatment of canopy complexity will be crucial to model performance in these aging forests. In addition, the relatively crude state of C-N integration (Yang et al. 2009) in most current models means that foliar N reallocation prompted by disturbance (H1.3) will be difficult to replicate, particularly for the gap-based model. These models will thus vary in their fidelity to observations, implying that the subtle disturbances of aging forests have fundamentally different effects on forest function than the stand-replacing disturbances of most ecosystem models (e.g., Peng and Apps 1999, Running 2008).

Figure 3. Model skill in replicating FASET results. Three very different models exhibited 20-80% declines in forest production, in sharp contrast to observations (dark line). Bond-Lamberty, Gough et al. 2015, Biogeosciences.

O3: Using newly available observations and a suite of ancillary data, broadly characterize age-NPP and NEP trajectories separately for temperate deciduous and coniferous forests.

  • H3.1: An empirical synthesis will reveal shallower post-peak NPP and thus smaller NEP declines in North America’s temperate deciduous forests than in coniferous forests, as well as interactive effects between disturbances (severity and frequency) and leaf habit.

    Luyssaert et al. (2008), in a transformative but controversial paper nearly a decade ago, reported NEP declines of 50% from peak through 100 years using an aggregated dataset of temperate and boreal forests, while postulating ecosystem structure’s importance in supporting production in older forests. We expect North American temperate deciduous forests, which are increasingly subjected to moderate rather than severe disturbance (Cohen et al. 2016), will exhibit shallower post-peak NEP declines. Indeed, an Ecosphere Innovative Viewpoint authored by the PIs (Gough, et al., 2016) suggests lesser NEP declines for deciduous forests of 16% through the century mark (Fig. 4). While deciduous and coniferous forest production may respond similarly to a given level of disturbance intensity (Hudiburg et al. 2009, Hardiman et al. 2013b), we expect that regional differences in disturbance regime, in part, shape age-NPP-NEP trajectories (Williams et al. 2014, Cohen et al. 2016), with North American deciduous and coniferous forests generally occupying very different space along the disturbance severity continuum (Gough et al. 2016). No synthesis explicitly contrasting temperate forest types has been conducted and newly available (e.g. FLUXNET) site-level flux (NPP and NEP), disturbance, climate, and ecological data make this analysis timely and critical; ~90% of extant data has been released since Luyssaert et al. (2008).

Figure 4. Our synthesis of temperate deciduous forests only suggests shallower NEP declines than expected from syntheses dominated by coniferous evergreen forests. Gough, Bond-Lamberty et al. 2016, Ecosphere.

This work’s open data and reproducible code (cf. our soil respiration database, Bond-Lamberty and Thomson 2010, that has been used and/or cited by ~150 studies, as well as the current “open experiment” at https://github.com/bpbond/cpcrw_incubation), will be a critical step and data platform for advancing a key unresolved area of ecosystem science: how and why ecosystems differ in their production trajectories (Ryan et al. 2010, Drake et al. 2011). O3 thus ties directly to broader impacts below, supporting sciencebased management tools to foresters (3) and providing resources for new analyses, beyond the immediate project and by other researchers (4).

III. Proposed Research and Methods

By tightly coupling field, model, and data experiments, our three-pronged approach aims to elucidate biological controls on functional resilience at both the ecosystem and regional scales, understanding why models have trouble simulating the consequences of moderate forest disturbance.

a. Site-level Manipulations of Disturbance Severity

The experimental portion of our proposed work will be conducted at the University of Michigan Biological Station (UMBS, 45 35.5’ N, 84 43’ W), in a secondary successional forest that is comprised of bigtooth aspen (Populus grandidentata), northern red oak (Quercus rubra), red maple (Acer rubrum), paper birch (Betula papyrifera), and eastern white pine (Pinus strobus). Average overstory tree age is ~100 years, decades older than the majority of forests in the U.S. northeast and upper Midwest. The plant community composition is typical of many forests in the upper Great Lakes region (USDA, 2008).

A principal research objective of NSF and DOE-supported long-term UMBS forest ecosystem studies, on which Gough is co-PI and Bond-Lamberty is a collaborator, is to identify physical and biological mechanisms controlling C cycling in disturbed and aging mixed deciduous-conifer forests. The Forest Accelerated Succession ExperimenT (FASET), in which >6,700 Populus and Betula trees were stem girdled in 2008 within a 39 ha area, has investigated how C pools and fluxes change as forests undergo a transition from early to later successional canopy dominants. The overarching FASET hypothesis, supported by results from this landscape manipulation, is that forest production across much of the region will be resilient following partial canopy defoliation and subsequently increase as canopies become more biologically and structurally complex, and as nitrogen (N) not taken up by senescing aspen and birch trees is redistributed to other species assuming canopy dominance. The experiment employs paired C cycling measurements within separate treatment and control meteorological flux tower footprints. The girdling treatment expedited mortality of early successional aspen and birch, promoting an emerging canopy that structurally and compositionally approximates projected ~50-yr regional successional changes (e.g., Wolter and White 2002). FASET’s experimental disturbance is similar in severity to a single level of moderate disturbance from wind, insects, and pathogens (Amiro et al. 2010, Hicke et al. 2012). Support from a DOE Ameriflux Core Site award ensures continuous collection of C cycling and meteorological data from flux tower footprints during the duration of work proposed here.

FASET was a crucial first experiment in advancing ecological understanding of the C cycle’s response to moderate severity disturbance, but it was designed to examine landscape scale C cycling processes following a single level of disturbance severity. Inherent spatial variation in tree mortality, a function of aspen and birch density, allowed for intriguing characterization of wood NPP’s response to multiple levels of disturbance severity (Fig. 1), but the experiment did not systematically characterize mechanisms underlying forest production’s resilience across a range of disturbance severities, or link directly with planned modeling experiments (Nave et al. 2011, Gough et al. 2013). By contrast, the work proposed here tightly integrates empirical and modeling experiments throughout, and explicitly links stand-level mechanisms to larger-scale analyses. We stand on the shoulders of FASET, in particular its successful approach to manipulating and quantifying the C cycle’s response to disturbance, but aim to unambiguously identify the mechanisms driving forest resilience or decline, and understand how and why models may fail to replicate observations.

b. Experimental Layout

As a new, integral component of this landscape-scale study, we propose a replicated ecosystem manipulation of disturbance severity to systematically identify thresholds and mechanisms associated with either NPP resilience/resistance or decline (Fig. 5). We will implement disturbance severity treatments that target LAI declines of 45 % (~FASET’s peak defoliation level), 65 % and 85 % in addition to an unmanipulated control. This range of disturbance severity encompasses and significantly extends FASET’s disturbance level, permitting us to identify the level of disturbance and associated mechanisms that cause NPP to decline. We will install four replicates of the four disturbance levels, for a total of 16, 0.5 ha treatment plots, each containing approximately 450 canopy trees (DBH ≥ 8 cm). Using 10-15 yr UMBS inventories and assessments of soils and plant communities, we have identified four replicates within a middle range of forest productivity (NPP of ~6 Mg C ha-1 yr-1). We limit our plots to areas outside of the C flux footprint of a meteorological (eddy-flux) tower in the unmanipulated, control forest (Gil Bohrer, UMBS meteorologist, personal communication).

Figure 5. Proposed experimental layout, with 4 full replicates of disturbance severity (expressed as fraction of canopy defoliation). The large semicircular plot and 3, 2 ha square plots (gray pixelated) illustrate locations of FASET plots installed in 2008 to examine landscape level responses to moderate disturbance that caused a temporary 44 % decline in LAI.

Each 2 ha treatment replicate is a grouping of four, 0.5 ha circular plots, each with a 25 m unmanipulated forest buffer, in which our defoliation treatments were randomly assigned. Circular plots comprising a replicate are irregularly clustered to control for inherent spatial variation in productivity (i.e., NPP), soils, tree communities, topography, and to avoid property boundaries and existing long-term experimental plots. Two random 0.1 ha subplots (together containing ~180 canopy trees) will be established in each plot interior for C and N cycling, physiology, phenology, and microclimate measurements (below). The UMBS Land Use Committee and Director (K.J. Nadelhoffer), and FASET PI (P.S. Curtis) have approved the location of proposed treatment plots shown in Fig. 5.

The main experiment (years 2-4) will follow an initial field season of pre-treatment data collection (year 1). We will systematically impose a gradient of disturbance severity by stem girdling dominant canopy trees, chosen to meet defoliations of 45 to 85 % (in addition to the control). We will use site- or region-specific allometry relating stem diameter to LAI to meet targeted defoliation levels. Early-successional aspen and birch will be marked for girdling first, with additional canopy dominant trees (of other species) randomly selected for girdling thereafter to reach targeted disturbance severities. During April-May 2017 (when the phloem and xylem are easily separated) and following one full year of pre-treatment data collection, we will employ sawyers (as in the successful FASET experiment) to score targeted trees and use pry bars to remove bark. This will simulate a range of disturbance intensities commonly prompted by insect outbreak, wind, or disease (Frelich 1995, Frelich and Reich 1995, Chen and Poland 2009, Nuckolls et al. 2009, Ford et al. 2012, Cohen et al. 2016). Stem girdling effectively killed aspen and birch trees in the first phase of the FASET experiment (Nave et al. 2011), resulting in rapid physiological decline, on average, within 1 year and mortality within 2-3 years. Other studies demonstrate a similarly rapid, often 1-year, return to pre-disturbance NPP following moderate disturbance (Amiro et al. 2010). We thus are confident that our 3 years of post-treatment measurements will coincide with peak LAI losses and, importantly, the critical period of dynamic canopy structural reorganization that we hypothesize constrains the trajectory of C cycling resilience or decline (see Fig. 1). Girdling eliminates transport of photosynthate to roots via the phloem, but allows xylem water transport. Unlike trunk removal, root suckering is minimal because root-stimulating auxins do not accumulate.

c. Research Methods – Interaction between Field and Modeling

Sections (d) and (e) detail the field and modeling methods, respectively. These activities will be closely coupled through a feedback loop between sampling and modeling (Fig. 6), which takes inspiration from recent calls to improve model-experimental communication and productivity (Dietze et al. 2013, Medlyn et al. 2015). Our aim is to jointly strengthen fundamental empirical and modeling knowledge of the C cycle’s response to moderate disturbance. The experimental design of the disturbance gradient (above), the measurements and their sample sizes (below), and the initial model experiments (below) are based on our biological understanding of this ecosystem, the sampling requirements to quantify NPP robustly (Gough et al. 2008), and our a priori knowledge of the models’ strengths and weaknesses. In each year, model assimilation experiments will iteratively inform the next field season’s sampling priorities, in particular by adjusting sample sizes. Further details are provided in the relevant sections below.

Figure 6. Overall project design, in which the initial field measurements are informed by a priori understanding of UMBS biology from the FASET project and model characteristics. After each field season, we interrogate the models (and potentially results of data synthesis): what adjustments will be most valuable to make the following year?

d. Research Methods – Field [Hypotheses addressed in brackets] We will conduct C and N cycling measurements to quantify total NPP and identify the mechanisms that enable resistance to moderate disturbance. Inventory approaches across a full range of disturbance severities will be used to estimate NPP in each of the two subplots nested within the 16 treatment plots. Even with disturbance-related increases in forest structural heterogeneity, NPP estimates at UMBS are made with high certainty following disturbance (95 % C.I. of estimated total NPP averages ± 16 %) (Gough et al. 2007a, Gough et al. 2010, Gough et al. 2013), and closely parallel independent meteorological estimates of NEP (Gough et al. 2013). We will calculate wood NPP as the annual sum of total subplot wood mass increment per unit area, determined from band dendrometers fitted permanently to a quarter of all live (prior to treatment) trees ≥ 8 cm stem diameter (D, n = ~700 trees total from all plots). Pre-treatment D will be recorded for all trees within each subplot and the D increment of trees without dendrometers inferred from fitted trees of the same species and disturbance level (Gough et al. 2010). We will corroborate inferred estimates by re-measuring the diameters of all trees in year 3, including those that have newly grown into ≥ 8 cm D class. We will estimate wood mass using site- and region-specific allometric equations relating D to wood mass (Gough et al. 2008). Leaf NPP will be estimated from annual litter mass production. Annual overstory leaf production will be quantified using five, 0.25 cm2 litter traps per subplot. Understory production will be estimated from two, 4 m2 plant survey quadrats within each subplot to determine fern and seedling/sapling production of individuals with D < 8 cm (Gough et al. 2007a). We will quantify fine root NPP from 5 root in-growth cores (1 m depth) harvested twice annually from each subplot, a frequency and sample size that UMBS minirhizotron data show is sufficient for accurately quantifying fine root growth and turnover before and after disturbance (Gough et al. 2008, Nave et al. 2011). Components of NPP will be summed to derive total NPP and error estimated according to Gough et al. (2008). [H1.1, 2.1]

We will quantify foliar N mass and reallocation among senescent and undisturbed residual plants, including those in the understory, and estimate whole canopy N mass via annual sampling of green-leaf N and leaf litter mass (Dronova et al. 2011, Nave et al. 2011, Gough et al. 2013). We will estimate foliar N reallocation and retention, rather than quantify N losses directly, because canopy N at our site is closely coupled with physiological properties that constrain NPP (Goodrich-Stuart et al. 2015). Canopy N losses from leaching, denitrification, and microbial immobilization, even in plots experiencing high levels of disturbance, were physiologically negligible within FASET’s moderate range of disturbance (Nave et al. 2011, Gough et al. 2013, Nave et al. 2014), but may be a mechanism causing production to decline following more severe disturbance (White et al. 2004, Gough et al. 2007a); therefore, we will carefully track and infer canopy N losses from inter-annual declines, should they occur, in total canopy N mass (Arain et al. 2006, Wang et al. 2013). Within one subplot in each of the 16 treatment plots, we will profile the green leaf N chemistry of the (~3 to 5) dominant species. Within each subplot, we will collect 20 leaves each from upper, middle, and lower (including understory) canopy positions annually during the peak of the growing season (July or early August). Dried leaves will be assayed for percent N on a Costech Analytical CHN analyzer (Costech International) in the UMBS Chemistry Lab. We will calculate species and whole canopy N stocks as the product of annual litterfall mass, determined from litter traps, and species’ foliar percent N values (Hardiman et al. 2013b). [H1.1, 1.3, 2.1]

We will track canopy phenology with a LAI-2200 Plant Canopy Analyzer (LI-COR Inc.) to determine whether the timing of leaf emergence and drop is altered by moderate disturbance (Potter et al. 2003). Readings will be taken every 3 m along 2 transects in each treatment plot on ~8 sampling dates from May (leaf emergence) to November (leaf drop). Phenology is a strong constraint on model performance, and a metric on which most models perform poorly (Richardson et al. 2012). In addition, a suite of leaf-level measurements will be used to quantify changes in physiological and morphological parameters of dominant and subdominant canopy trees and of “released” trees in the understory following disturbance (Gough et al. 2004), in particular maximum stomatal conductance (via LI-COR 6400), specific leaf area, and leaf C:N, which are highly influential parameters in most models including Biome-BGC and ED2 (e.g., White et al. 2000). [H1.2, 2.1]

We will quantify canopy structure, including metrics such as rugosity, gap fraction, and clumping indexes, annually along 2 perpendicular transects in each 0.5 ha plot at peak LAI using a portable canopy LIDAR (PCL). The design, operation, and validation of this system is described in Parker and Russ (2004) and was used previously at our site, following the development of analytical software, to relate canopy structural complexity to NPP (Hardiman et al. 2011, Gough et al. 2013, Hardiman et al. 2013a, Hardiman et al. 2013b) (Fig. 2). The PCL is based on an upward looking, near-infrared pulsed-laser operating at 2000 Hz (model LD90-3100VHS-FLP, Riegl USA). Prior analysis of UMBS forests indicates that canopy structural metrics derived from LIDAR measurements may be stronger correlates of stand-scale productivity than LAI; they also provide precise empirical height profiles, a much more stringent diagnostic test for model performance than simply LAI. [H1.1, 1.3, 2.1] We will quantify disturbance-related shifts in canopy light distribution, and soil temperature and moisture using an array of solar and battery powered sensors deployed in one subplot within each of the 16 treatment plots. PPFD will be recorded to a datalogger at 0.1 Hz and averaged over 10 minute intervals during the leaf-on period with an array of eight randomly distributed sensors (LI-190 Quantum Sensor, LI-COR Inc.), and used to determine growing season fAPAR and APAR (Beaudet et al. 2007, Kranabetter and Simard 2008, Tobin and Reich 2009). An existing BF5 sunshine sensor (Delta-T Devices, Cambridge, UK) mounted on a nearby tower measures above canopy total and diffuse radiation for the estimation of fAPAR. Year-round soil temperature (7.5 cm) and moisture (20 cm) will be recorded hourly using a datalogger in 5 randomly assigned locations within each of the 16 plots via Hobo temperature (thermocouple) and moisture (TDR) probes (Onset Computer Corporation, Pocasset, MA USA). [H1.1, 1.2, 2.1]

e. Research Methods – Modeling

The modeling component of the project will parallel and complement the fieldwork, asking when NPP declines with increasing disturbance, why the models fail (as we hypothesize) to match field observations of NPP and nitrogen reallocation, and, as a byproduct of these questions, which field-measureable parameters are critical to constrain these model components. We will use two complementary models running under a unified data assimilation framework. The first is a version of Biome-BGC, a classic “big leaf” model with coupled water, carbon, and nitrogen cycles (Running and Gower 1991, Thornton and Rosenbloom 2005), and a Farquhar photosynthesis submodel. Its algorithms form the basis of NCAR’s Community Land Model (CLM), and the model has seen significant recent development work (Higy et al. 2016). The second is ED2, a model that uses size- and age-structure partial differential equations (PDEs; Moorcroft et al. 2001) to approximate the behavior of a stochastic gap model at medium to large scales. ED2 has been used for a variety of optimization and data assimilation exercises (Medvigy et al. 2009), and was recently rewritten by Dietze and Matthes (2014) to investigate the impact of FASET-like disturbances. The complementary nature of these models is critical: one is a “big leaf” biogeochemical model focusing on process representation, while the other emphasizes mathematical scaling of a gap model across time and space; both are part of the development stream of CLM, as well as DOE’s new Accelerated Climate Modeling for Energy (ACME). This provides a strong framework for examining if experimental treatment effects can be reproduced in diverse modeling experiments.

Model and data development. We will use the NSF-funded PEcAn framework for ecosystem model analysis and data assimilation (http://pecanproject.github.io). PEcAn is an open-source, modular workflow that manages the flow of information into and out of ecosystem models, providing a set of tools to quantify and propagate parameter uncertainty (e.g., Dietze et al. 2014). Prior probability distributions for model parameters are stored in the PEcAn database and automatically updated with a hierarchical Bayes meta-analysis when new data are added to the database. These updated distributions are sampled to drive an ensemble of model runs, which can quantify parameter sensitivity and model prediction uncertainty. PEcAn currently supports ED2, and given our experience with Biome-BGC (Bond-Lamberty et al. 2005) and PEcAn support (see letter of support), we do not anticipate problems adding this model into the PEcAn framework. PEcAn’s primary PI, Dr. Michael Dietze, has indicated that single-site version of CLM 4.5 (Oleson et al. 2013) will be added during our project’s funding period, and we will test it as well. Biome-BGC already has an existing disturbance mechanism that can, with minimal effort, be modified to include the effects of stem girdling on canopy structure as proposed here. ED2 models subgrid (~10 ha) disturbance heterogeneity using its PDEs to approximate the behavior of a spatially distributed ensemble of individual plants.

Model initialization and optimization. The first phase of our modeling work will optimize the models by using (i) the known parameter sensitivities of the two models, (ii) eight years of post-treatment FASET data, to which Gough has full access as project co-PI, and (iii) ongoing measurements in the unmanipulated control site (with the selection of UMBS as a long-term Ameriflux Core site ensuring C flux data from treatment and control sites through at least 2020). The breadth and length of these records provide a strong constraint for the modeling experiments, as meteorological variability means that data assimilation experiments greatly benefit from longer records (Hill et al. 2011). Model optimization is a standard procedure using PEcAn (Dietze et al. 2014). Critical parameters–e.g., C:N ratios, specific leaf area–have been well characterized at UMBS. In summary, this phase will result in model parameter sets that produce the best fit between model outputs and decadal (pre-experiment) observations at UMBS.

Model testing. The second modeling phase will examine at what point, and why, the models fail in replicating how NPP responds to the three disturbance intensity treatments. By the end of year 3 (2018, see timeline below) we will have two years of data on forest response to experimental treatments against which the models will be tested. Following the hierarchical approach of Wang et al. (2009) and the recommendations of Medlyn et al. (2015), we will use three functionally cascading tiers (leaf level, plant-level C and N cycling, and litter/soil dynamics) for analysis. These tiers are mostly common to both models (except for N cycling, which will be Biome-BGC only) and map directly onto the hypothesized structural and allocation reasons for ecosystem resilience. Within each tier, we will quantify model deviation from observations based on root-mean-square difference between modeled and observed data (Taylor 2001), to identify the model tier and then specific processes most responsible for the models’ hypothesized failure to simulate NPP resilience, increased structural complexity, and N reallocation. The models will be run across the full suite of disturbance severity within the PEcAn framework. [H2.1] Uncertainty estimation. We will focus on uncertainty estimates from two sources: the observational means and errors associated with critical parameters, and interannual climate variability effects from ensembling (Thornton et al. 2002). We will use these results to estimate parameter means, distributions, and correlations, providing insight into model function, measurement priorities, and improved predictive understanding. At the moment we are particularly interested in using PEcAn, Biome-BGC and ED2 to “rate the data” (Keenan et al. 2013a), to determine which measured variables are relatively redundant and which should be prioritized, i.e., be given additional resources (increasing sample sizes) the following year. We will perform this exercise annually (Fig. 6). [H2.1]

f. Research Methods – Synthesis of NPP and NEP-Age Trajectories in Temperate Forests

Following O3.1, the synthesis will use previously unavailable data to examine how temperate forest C cycling changes with age, disturbance frequency and severity, and leaf habit in a fully reproducible framework. We estimate, based on our recently published synthesis of deciduous forests only (Gough et al. 2016), that ~90% of extant data have been released since the well-known analyses by Luyssaert et al. (2007, age-focused) and Amiro et al. (2010, disturbance-focused), making our analysis of interacting age-disturbance effects timely and needed. We focus on North America’s temperate deciduous and coniferous forests since past and present disturbance regimes are well differentiated for these forest types (Cohen et al. 2016, Gough et al. 2016), providing a rigorous test of H3.1.

We will perform two complementary syntheses. The first uses a “space-for-time” substitution employed in several highly influential age-production analyses (Luyssaert et al. 2008, Tang et al. 2014), and the second a site-level time-series analysis (Amiro et al. 2010) that explicitly couples interannual and decadal changes in NPP and NEP with disturbance. The first approach allows examination of broad multi-century trends in production as shaped by disturbance, but is limited by potential co-varying drivers of production and similar artifacts common to multi-site analyses (cf. Magliocca et al. 2015). The second, site-level, analysis takes advantage of a growing number of mid-successional temperate forests with 10+ years of continuous NPP and NEP data (N = 7 deciduous, N = 4 conifer based on the current FLUXNET2015 database) to examine whether disturbance modifies the post-peak production trajectories of a collection of mid-successional forested sites, similar in approach to recent analyses examining environmental drivers of decadal production (Keenan et al. 2013b). This will provide new explanatory power for interpreting trends reported by Amiro and co-authors: variable, but unattributed, shifts in NEP among forest ecosystems subjected to moderate disturbance (see Amiro et al. 2010, figure 5).

Both space-for-time and time-series syntheses will use recently released (Dec. 2015, July 2016) FLUXNET2015 Tier 1 (open and free use) and Tier 2 (collaboration/consultation required) NPP and NEP data (the latter provided in Biological, Ancillary, Disturbance and Metadata, BADM) and, because not all site data are posted via FLUXNET, published NPP and NEP estimates for North American temperate forests. Presently, 18 and 57 North American temperate deciduous and coniferous forests, respectively, are included in the FLUXNET2015 release (http://fluxnet.fluxdata.org/data/data-availability/), with upcoming releases including additional sites. Also, we have assembled data for 65 North American temperate forests (as part of Gough et al., 2016 synthesis); some, but not all, overlap with FLUXNET2015. The FLUXNET2015 release increases sample size, will ensure standardized processing of eddy covariance-based NEP, and provides key ancillary climate and disturbance data. We will supplement the FLUXNET2015 and Gough et al. datasets with ancillary data–in particular, the CRU (Climatic Research Unit, U. of East Anglia) historical climate data, the SoilsGrids1km three-dimensional soils database (Hengl et al. 2014), and ORNL DAAC N deposition data (Dentener et al. 2006).

Site disturbance – severity and frequency – will be derived from TimeSync (Cohen et al. 2010, Cohen et al. 2016), and compared with on-the-ground disturbance records when available. TimeSync is a highly calibrated and validated (Schroeder et al. 2014) visualization tool for quantifying disturbance severity at a high resolution via change detection of 30m x 30m Landsat time-series imagery. The coherent timespan and fine spatial scale of TimeSync’s inferred disturbance data, coupled with FLUXNET’s NPP and NEP data, will allow for a more quantitative analysis of the effects of disturbance on C cycling than previous efforts (e.g. He et al. 2012). Several sites, including (deciduous) Harvard Forest, UMBS and (coniferous) Metolius, have published records of disturbance severity and frequency (Urbanski et al. 2007, Meigs et al. 2009, Gough et al. 2013). These observational data are important as a validation of TimeSync data, and because they show that disturbance-related tree mortality (e.g., comparable to NPP in some years; Urbanski et al. 2007) and LAI decline (e.g., 44 % reduction from peak; Gough et al. 2013) is significant for many sites and well within the disturbance severity and spatial scale detection limits of Landsat and TimeSync (Cohen et al. 2010, Kennedy et al. 2012, Pflugmacher et al. 2012).

Because the FLUXNET data are not amenable to formal meta-analysis techniques (Hedges et al. 1999), we will instead use two broad regression and machine-learning tools in both the space-for-time and time-series approaches. The first uses non-parametric and frequentist change-point detection algorithms (Killick et al. 2012), implemented in the changepoint and BreakoutDetection packages in R, to identify potential transition phases during succession for forest stands with similar disturbance characteristics (in space-for-time analysis) and changes in site production trajectories following disturbance (in time-series analysis). These algorithms are entirely data-driven (i.e., they do not require a pre-specified number of segments or inflection points). We will also compute standard metrics such as quantity and time of maximum C fluxes. In a second analysis, we will examine how different factors affect these change-points, transition phases, and peaks using Random Forest (RF), a nonparametric machine learning technique for classification and regression. RF calculates variable importance by aggregating regression trees using different random samples of the data (Breiman 2001). This will allow us to quantify how temperate forest NPP and NEP trajectories differ with disturbance, region, climate, soils, N deposition, and leaf habit. The PIs have extensive experience with these techniques (Bond-Lamberty et al. 2012, Bond-Lamberty et al. 2014, Stuart-Haentjens et al. 2015).

This work will produce syntheses, as described above, but also an open platform for future analyses. The PIs have a long and successful history of collaborative and synthetic analyses (Bond-Lamberty et al. 2015, Gough et al. 2016) that have generated nearly 1000 citations and catalyzed additional research and workshops. The PIs have also pushed strongly for open data and reproducible science (Rueegg et al. 2014, Bond-Lamberty et al. 2016): PI Gough contributes to the FLUXNET datasets described above, while Bond-Lamberty has contributed to, and curates, several databases (Bond-Lamberty and Thomson 2010, Kim et al. 2012, Falster et al. 2015) that many other scientists have used and built upon for subsequent analyses. This experience gives us the motivation and skills to create an open “platform” – a code repository for reproducible data processing (that can be easily re-run after new FLUXNET data releases, for example); web-accessible data visualization tools, using the shiny framework (http://shiny.rstudio.com); and reproducible and transparent analyses and manuscripts (the latter using RMarkdown). This will extend and improve on efforts, by both FLUXNET and the broader science community, to understand age-production successional trajectories, and serve as a tool and resource readily accessible to researchers, educators, and policy-makers.

g. Testing the Hypotheses

  • H1.1: We will quantify annual nitrogen- and light-use efficiencies as NPP per unit canopy N mass or light absorbed (i.e., APAR), respectively (Runyon et al. 1994, Gower et al. 1999, Dronova et al. 2011, Hardiman et al. 2013b). We will quantify annual peak LAI, and conduct annual LIDAR scans of the canopy to derive quantitative metrics of canopy structure (Hardiman et al. 2011, 2013). We will examine whether shifts in resource-use efficiency over time correspond with disturbance-generated changes in canopy structure, and determine the threshold at which greater resource-use efficiency no longer compensates for LAI losses.

  • H1.2: We will compare light distribution between, and leaf-level physiology and morphology of, the canopy dominant and subdominant vegetation as structural shifts unfold following disturbance. We expect that leaf-level photosynthesis of subdominants will increase significantly as canopy gaps form and light penetrates deeper into the canopy, while photosynthesis of dominant canopy individuals remains relatively stable; this will drive an increase in ecosystem light-use efficiency. We will initially focus on SLA and leaf N, critical parameters for scaling and ecosystem modeling (Reich et al. 2007), but these may change depending on the data assimilation experiment results (described above).

  • H1.3: We will examine changes over time in species’ foliar and whole-canopy N mass, and assess whether complete N reallocation in support of LAI replacement together with enhanced canopy structural complexity (e.g., higher gap fraction and rugosity, lower clumping index) confer greatest (i.e., most rapid and complete) resilience or resistance to moderate disturbance.

  • H2.1: After optimizing the parameters to which each model is most sensitive, we will evaluate model skill, propagating both observational (plot-to-plot) and model-derived (as described above) errors. This will be performed in the PEcAn framework using several primary data streams such as NPP, leaf-level photosynthesis, canopy vertical profiles, and C:N ratios. For example, we will test how the models’ fixed C:N ratios constrain their performance by quantifying how close the models’ optimum canopy N parameters are to observed values. We anticipate they will diverge as the imposed disturbance intensity increases, and that ED2 in particular will suffer from its lack of a nitrogen cycle.

  • H3.1: We will test our a priori hypothesis, that deciduous forests will exhibit shallower post-peak NPP and NEP declines than coniferous forests, against the results of the data-driven synthesis. The algorithm-derived variable importance rankings will be tested for spatial variability and influential outliers, and a pseudo (out-of-bag) R2 calculated. Hypothesized interactive effects between disturbance severity, frequency, and leaf habit will also be tested. Finally, we will compare these synthesis results against the site-level model sensitivities derived from the data assimilation exercises, and assess to what degree UMBS is representative, or not, in the larger landscape of similar forests.

IV. Broader Impacts of the Proposed Work: Integration of Research, Education, and Outreach

In addition to the broader impacts detailed below, the work will advance a productive, but to date unfunded, collaboration between middle career scientists at the Joint Global Change Research Institute and Virginia Commonwealth University. Additionally, this project will contribute to a prolific legacy of forest ecology experiments at UMBS, and support present and future research and teaching at the Station.

1) Enhance undergraduate research training opportunities for VCU’s diverse student body.

We propose a new undergraduate research training program, with an emphasis on underrepresented students in the sciences, which will support student research activities at UMBS and academic year research mentoring at VCU. The benefits of such training to students of the sciences are well-established, and include an increased capacity to formulate hypotheses and execute research activities, analyze data, and effectively communicate research results (Kardash 2000, Landrum and Nelsen 2002, Bauer and Bennett 2003, Seymour et al. 2004, Hunter et al. 2007). Undergraduate students with graduate preparatory experience are also more likely to pursue Ph.D. degrees in the sciences (Hathaway et al. 2002, Russell et al. 2007, Junge et al. 2010). This trend is particularly true among underrepresented groups in the sciences, with effective undergraduate research training programs often providing transformative career-shaping experiences (Barlow and Villarejo 2004, Villarejo et al. 2008, Hurtado et al. 2009, Jones et al. 2010).

VCU enrolls and awards a high percentage and volume of baccalaureate biology degrees to ethnic and racial minorities, making the program an ideal target for enhancing the graduate preparedness of underrepresented minorities in the sciences. VCU Biology trains the largest number of undergraduate majors at the university (1,974 of 22,335 students enrolled, 2014-2015), and is among the top 20 universities in the nation for biology degrees awarded to ethnic and racial minority students. In Spring 2015, 19.1 % of enrolled undergraduate biology students were African American, exceeding the national average of 13 % (NSF, 2013). A high percentage of biology undergraduates (26.6 %) were Asian/Pacific Islanders, a value that far surpasses the national average of 5.6 %.

We will use NSF support to cultivate and enrich ongoing efforts at VCU to enhance undergraduate research training by: a) initiating an annual “field research” day for undergraduates at VCU’s nascent field station, the Rice Rivers Center (http://ricerivers.vcu.edu/); b) directly engaging a subset of students in collaborative research at UMBS; and c) using data collected from UMBS to support a research training biology capstone seminar. We describe each of these:

  • a) Initiation of a “field research” day. Surveys of VCU biology undergraduates in an upper level field course instructed by PI Gough (Forest Ecology, BIOL 422) indicate that 55 % of students enrolled are pre-health and 20 % are undecided in their career choice (n=48). These surveys additionally suggest that students would consider non-professional science careers, though many are unaware of other options. We will use NSF support to initiate a “field research” day at VCU’s field teaching and research station (located 50 km east of main campus), in which 40 undergraduate first-year students enrolled in an existing Introduction to Biological Science I (BIOL 151) course will participate in an all-day field research workshop. The intent of the workshop is to inform and educate biology students in their first year of the breadth of ecological research opportunities available to prospective undergraduate researchers and to showcase career options in ecology. The workshop will highlight aquatic, terrestrial, and wetland research conducted at the Rice Rivers Center. The PIs will open the workshop with an ecology career seminar. Post-workshop surveys will be conducted to determine the effectiveness of the workshop in enhancing understanding of career options in ecology.

  • b) Provide summer fellowships for four VCU undergraduates to attend UMBS as research investigators. Student collaborators will reside at UMBS during the 8-wk summer session and serve on the broader Forest Ecosystem Study team, which perennially includes 20-30 students, postdocs, and PIs from multiple institutions. Our research group has trained scores of undergraduate investigators through the long-term NSF supported REU program titled “Biosphere-Atmosphere Interactions in a Changing Global Environment”. PI Gough has mentored 17 REU participants, with 10 undergraduate advisees co-authoring 6 published manuscripts (Gough et al. 2007a, Gough et al. 2007b, Gough and Elliott 2012, Hardiman et al. 2013b, Fahey et al. 2016, Schmid et al. 2016); one manuscript with 2 undergraduate co/authors (one as lead) is in review. Given VCU’s diverse student body, we anticipate that directed recruiting of underrepresented groups will be unnecessary. However, we will advertise through the VCU Undergraduate Research Office, and model our fellowship application after that for the UMBS NSF REU program, which also seeks to enroll underrepresented students in the sciences.

  • c) Data in support of a biology capstone seminar. Undergraduate students interested in pursuing advanced degrees in the sciences are often ill-prepared to navigate the graduate application and admissions process (Coronado et al. 2010). UMBS data collected by VCU undergraduate research investigators will be summarized and presented to peers during a new biology capstone seminar intended to prepare graduating seniors for the (non-professional) graduate school application process. VCU students who have conducted undergraduate research and who are pursuing graduate training will be targeted. The course, directed by PI Gough, will require undergraduate student researchers to present their findings (including participants in UMBS research), assist students in identifying a suitable graduate program in the sciences, and will offer approaches to applying to graduate school. Course topics will include how to appropriately contact candidate graduate mentors, seek student financial support for graduate school, and effectively prepare for the Graduate Record Exam (GRE).

2) Pursue educational and collaborative tie-ins for high school teachers and students. We will work with primary and secondary educational institutions, taking advantage of the DOE Science Undergraduate Laboratory Internship (SULI) at PI Bond-Lamberty’s institution, and Workforce Development for Teachers and Scientists (WDTS) programs. The PIs, who have existing roles in K-12 science education, will work with educators to develop forest C Cycling laboratory modules for elementary through high school science classes.

3) Delivery of science-based management tools to regional foresters. We will host, at UMBS, a 1.5 day workshop for natural resource managers during the autumn of the fourth funding year titled “Complexity in Forests of the Upper Great Lakes: Merging Science and Management”. The objective of the workshop will be to deliver to foresters of the upper Great Lakes region science-based strategies for enhancing ecosystem goods and services, including forest C storage, by implementing management practices that augment forest structural complexity. This workshop will build on emerging science-driven silvicultural paradigms that aim to diversify services derived from forests by increasing, rather than diminishing, forest structural complexity (Puettmann et al. 2009). UMBS and its research faculty, including PI Gough, have a strong record of hosting and delivering science-based management tools to foresters. Gough maintains professional relationships with natural resource-focused NGO, USDA Forest Service, and Michigan DNR personnel, and will advertise the workshop broadly among these groups. He is also a past member of the Society of American Foresters (SAF) and will work with this organization to ensure that participating foresters receive credit in support of their SAF Certified Forester certificate. During the workshop, we will offer tours of the experimental site, conduct breakout sessions to identify common science and silvilcutural knowledge gaps in relation to the workshop’s theme, and deliver science-based, including modeling, management tools through presentations, technical documentation, and demonstration modeling exercises. Biome-BGC has seen significant application in managed forests (e.g., Cienciala and Tatarinov 2006), and we will use an Excel version of it for hands-on demonstrations.

4) Data and code platform for future age-NPP and NEP analyses. The synthesis component of the project (section IVf), leverages our extensive experience with reproducible methods and open data (e.g. Bond-Lamberty and Thomson 2010, Bond-Lamberty et al. 2016) to lay the foundation for future modeling and empirical analyses addressing how and why age-NPP/NEP trajectories differ among forest types. Though our emphasis is on ecosystem structure, we expect such a resource (of fully reproducible code and data) to be broadly useful, providing significant “value added” to FLUXNET and other data sources.

V. Personnel

Gough and Bond-Lamberty have a long-standing working relationship, recently co-authoring papers that examined modeled responses to moderate disturbance (Bond-Lamberty et al. 2015) and marshaled evidence for sustained production in aging temperate deciduous forests (Gough et al. 2016). The two postdoctoral researchers, one hosted by VCU and the other at PNNL, will be essential team members, participating in data collection and analysis, modeling activities, mentoring, and dissemination of our research results (see Postdoctoral Researcher Mentoring Plan). Gough will advise one VCU M.S. student annually on the project (2 over the course of the project), with Bond-Lamberty on the students’ committees, and Gough will additionally supervise one undergraduate summer research collaborator annually at UMBS and during the academic year. NSF’s responsible conduct of research (RCR) training will be fulfilled by all project-affiliated undergraduate and graduate students, and postdocs, through their successful course completion of Responsible Scientific Conduct (OVPR 602). PI biographies:

Christopher M. Gough, Virginia Commonwealth University, Research Assistant Professor: Gough will serve as project manager, co-direct field and modeling activities, and advise students and a postdoc. He is a forest ecologist, a 13-year UMBS investigator, and co/author of 55 peer-reviewed publications in journals including Global Change Biology, Ecology, Ecological Applications, BioScience, and Nature Education. Gough is also co-author of science policy documents, including the forthcoming Second State of the Carbon Cycle Report (SOCCR2).

Ben Bond-Lamberty, Battelle/JGCRI, Research Scientist: Bond-Lamberty will work with Gough on treatment implementation and measurements, leading a subset of field measurements as well as the modeling component of the project, and advise the PNNL postdoctoral scientist. He is a forest ecologist with extensive expertise in modeling ecosystem carbon cycling and synthesizing large empirical data sets (e.g., Bond-Lamberty and Thomson 2010) and has co/authored 64 peer-reviewed publications in journals including Global Change Biology, Oecologia, Ecosystems, Science, and Nature.

VI. Timeline

We request 4 years of support from NSF to conduct proposed research and outreach activities, with the first year (2017 committed to pre-treatment data collection. Broader impacts/outreach activities will occur throughout the funding period as specified below. Our proposed timeline is as follows:

Year 1 (2017): Treatment plot boundaries are marked and subplots located. Meteorological stations are installed in all 16 plots. Baseline C and N pool and flux, leaf-level physiology, canopy structure (i.e. LIDAR), and canopy phenology data are collected following installation of band dendrometers, in-growth cores, and litter traps. Simple modeling runs are performed with baseline data to check expected parameter sensitivities. Postdocs becomes familiar with the PEcAn toolbox and data assimilation (DA) system. Undergraduate research training program is initiated at VCU.

Year 2 (2018): Basic DA experiments incorporate year 1 measurements into Biome-BGC and ED2, quantifying model skill improvement and informing field sampling priorities. Targeted trees are girdled during spring, prior to leaf-out. Meteorological, C and N pool and flux, leaf-level physiology, and canopy structure and phenology data collected. C cycling module development begins for high school educators. Undergraduate research training program continues at VCU. FLUXNET and other ancillary data are acquired, and coding begins (in an openly available repository), in preparation for the data synthesis.

Year 3 (2019): DA incorporates year 2 measurements. Meteorological, C and N pool and flux, leaf-level physiology, and canopy structure and phenology data collected. Modeling experiments will test model sub-component performance against measured changes in canopy structure, N pool shifts, and major C fluxes. C cycling module finalized and distributed to high school educators. Undergraduate research training program continues at VCU. Synthesis and coding analysis proceeds, with sources updated as necessary when new FLUXNET and other data releases occur. Manuscript preparation emphasizes initial NPP/C cycling responses to moderate disturbance.

Year 4 (2020): Meteorological, C and N pool and flux, leaf-level physiology, and canopy structure and phenology data collected. DA experiments incorporate year 3 data and aim to understand what parameter sensitivities most critically determine performance against the final field data. Manuscript preparation continues, emphasizing the mechanisms causing NPP/C cycling resilience/resistance or decline following moderate disturbance, and the temperate forest synthesis work. Open synthesis data released and advertised (e.g., via ecolog-l, Amerflux listserves) following publication acceptance. Workshop for foresters is hosted by PIs at UMBS. VCU undergraduate research training program continues.

VII. Results from Prior NSF Support

Gough is co-PI (K. Nadelhoffer, PI) on an NSF LTREB award: “Drivers of forest C storage from canopy closure through successional time”. The award total is $448,585, with VCU receiving $80,543 over 5 years, 2014-2019. Following the second funding year, 3 publications are nearing completion or in review, including one authored by Gough’s MS student, and another co/authored by two NSF REU advisees. This project is relevant to work proposed here, examining how forest physical structure constrains decadal to century C cycling processes, from early through late stages of ecosystem development.

Bond-Lamberty is co-PI (with M. Harmon and R. Vargas) on a funded NSF Emerging Frontiers proposal: “Reducing Uncertainty in Heterotrophic Respiration: Linking Continental Experiments, Analytical Modeling, and Shared Data Sets”. The award total was $85,006, and funded two community workshops in 2013 and 2015. One resulting publication appeared in 2016, and one remains in progress.