In ashiklom/fortebaseline: Forest Resilience and Threshold Experiment (FoRTE) ED2 baseline runs source code and manuscript
library(knitr)
library(here)
library(magrittr)
library(kableExtra)
fig <- function(...) {
knitr::include_graphics(here::here("analysis", "figures", ...))
}
out_format <- knitr::opts_knit$get("rmarkdown.pandoc.to")
pgbreak <- function() {
if (grepl("docx", out_format)) {
return("##### pagebreak")
} else if (is_latex_output()) {
return("\\newpage")
}
}
opts_chunk$set(
collapse = TRUE,
warning = FALSE,
message = FALSE,
echo = FALSE,
comment = "#>",
out.width = "\\linewidth"
## fig.path = "../figures/"
)
## library(fortebaseline) # Or use devtools::load_all('.', quiet = T) if your code is in script files, rather than as functions in the `/R` diretory
n <- 0
mdcap <- function(text) {
n <<- n + 1
out <- sprintf("**Figure %d**: %s", n, text)
if (grepl("gfm", out_format, ignore.case = TRUE)) {
out
} else {
invisible(NULL)
}
}
m <- 0
htmlfig <- function() {
m <<- m + 1
out <- paste0("**Figure ", m, "**:")
if (knitr::is_html_output()) {
out
} else {
invisible(NULL)
}
}
Introduction
Human land use patterns over the last 500 years, and particularly over the last century, have led to decline in primary forests across most of the Earth's surface [@hurtt_2011_harmonization].
However, in many places, subsequent changes in land use (e.g. abandonment of agriculture and pasture) have caused much of this area to regrow as secondary forest [@hurtt_2011_harmonization].
Understanding the drivers behind this regrowth, and the ability to predict its future course, are therefore essential to shaping global policy around climate, biodiversity, and natural resource management.
The temperate forests of the Upper Midwest are a model system for studying secondary forest ecology and biogeochemistry.
These forests have been a net carbon sink for most of the past century, as fast-growing early successional species established and thrived following near-complete deforestation around the end of the 19th Century [@birdsey_2006_forest; @williams_2012_carbon].
More recently, many Upper Midwest forests are undergoing a shift in community composition from early-successional to mid-successional species, and the ecological and biogeochemical consequences of this transition are not well-established [@gough_2016_disturbance].
In addition, these forests face a multitude of direct and indirect pressures from both biotic and abiotic sources [@shifley_2016_future], including longer and more severe droughts [@gustafson_2016_implications; @swanston_2017_vulnerability; @lienard_2016_us] and non-native insects and pathogens [@lovett_2016_nonnative].
Collectively, these processes have culminated in forests that are different from their pre-settlement climax counterparts [@stearns_2002_one; @wolter_2002_recent; @thompson_2013_four], and our ability to make predictions about such non-analog environments based on the past is limited.
More reliable predictions can in theory be obtained by using dynamic vegetation models that explicitly represent processes involved in forest growth and mortality.
Vegetation models fall along a spectrum of complexity, with trade-offs between realism and process fidelity on one hand and computational demand and data requirements on the other.
On one side are "big leaf" models---such as PNET [@aber_1992_pnet], SiBCASA [@schaefer_2008_sibcasa], and Biome-BGC [@running_1993_generalization]---in which the vegetation at a particular location is represented by a single average "plant".
On the other are "individual-based" (or "gap") models---such as LANDIS [@mladenoff_2004_landis; @scheller_2007_design] and UVAFME [@shuman_2015_forest]---that explicitly simulate competition between multiple individuals [@shugart_2015_computer].
In the middle are models that approximate the emergent outcomes of competitive interactions among different individuals without explicitly simulating each individual [@purves_2008_predictive].
One example of this approach is the Ecosystem Demography model [ED; @moorcroft_2001_method; @medvigy_2011_predicting; @longo_2019_ed1], which is the focus of this study.
Regardless of vegetation model complexity, their projections are inherently uncertain.
This uncertainty comes from many different sources, of which we consider two in this study:
Structural uncertainty refers to which process are included and how they are represented.
Model intercomparison studies have shown the importance of different model structures to projections of overall land carbon sequestration [@friedlingstein_2006_climate; @friedlingstein_2014_uncertainties; @lovenduski_2017_reducing], response to CO~2~ fertilization [@zaehle_2014_evaluation; @medlyn_2015_using; @walker_2015_predicting], and soil C sequestration [@sulman_2018_multiple].
This structural uncertainty has been attributed to differences in model representations of key processes, including canopy radiative transfer [@fisher_2017_vegetation], soil biogeochemistry [@wieder_2017_carbon; @sulman_2018_multiple], stomatal conductance [@franks_2018_comparing], and photosynthesis [@rogers_2016_roadmap].
However, the extent to which specific processes contribute to model uncertainty is difficult to evaluate from model intercomparisons because models differ from each other in too many ways, meaning that structural effects are confounded with other uncertainty sources.
At the same time, parameter uncertainty arises because of natural variability and imperfect calibration in process representations.
Previous analyses of the ED model in particular have highlighted the importance of parameter uncertainty (and reduction thereof) to predictions of primary productivity and composition [@dietze_2014_quantitative; @raczka_2018_what; @viskari_2019_influence].
While many modeling studies have looked at structural and parameter uncertainty independently, comparative analysis of contributions of both uncertainty sources are rare in the ecosystem modeling literature.
Light availability and distribution drives demographic variation in plant C fixation and community structural changes through succession.
Observations show that competition for light drives niche differentiation [@silvertown_2004_plant] and consequently successional change [@canham_1990_light; @bazzaz_1993_successional], particularly in closed canopies.
Prior ED versions were limited in their capacity to sustain multiple plant functional types during stand development due in part to unrealistically aggressive shading of the understory and a lack of consideration for tree crown area, such that the marginal height advantage of a vegetation cohort resulted in a disproportionally large growth advantage [@fisher_2010_assessing].
Alternative approaches that accommodate horizontal competition for light [e.g. perfect plasticity assumption, @strigul_2008_scaling; @weng_2015_scaling] reduce the relative impact of height on light availability and growth, and may therefore promote competition [@fisher_2017_vegetation].
In addition, trait plasticity may permit individuals to acclimate and species to adapt to changing light environments in ways that improve light use efficiency [@niinemets_1998_analysis; @lloyd_2010_optimisation; @keenan_2016_global].
To this end, we focused our structural uncertainty analysis on three processes related to competition for- and utilization efficiency of light:
Representation of horizontal competition;
representation of canopy radiative transfer (especially of diffuse radiation);
and trait plasticity with light level.
This study focuses on the interaction between parameter and structural uncertainty.
Our study is organized around the following guiding questions:
(1) Which processes related to light utilization are most important to accurately model community succession and C cycling in Upper Midwest forests?
(2) What are the relative contributions of parameter uncertainty (data limitation) vs. structural uncertainty (theoretical limitation)?
We ran ensemble simulations of the Ecosystem Demography (ED-2.2) model with a factorial combination of submodels related to canopy radiative transfer formulation, horizontal competition, and trait plasticity.
We then used sensitivity and variance decomposition analyses to evaluate the contribution of parameter uncertainty for each model configuration, exploring the ecological implications of the results in the context of ongoing global change.
Methods
Site description
We performed this study at the University of Michigan Biological Station (UMBS; Ameriflux site US-UMd, 45.5625$^{\circ}$, -84.6975$^{\circ}$), located in Northern Lower Michigan, USA.
The area surrounding the research station is 87% well-drained upland forest and 13% wetland [@bergen_2007_observing]; the focus of this study is on the former.
The landscape geography of the UMBS upland forest is 20.4% moraine, 37.8% high outwash plain, 31.3% low outwash plain, 5.7% lake-margin terrace, 3.6% ice-contact, and 1.2% lowland glacial lake [@bergen_2007_observing].
The UMBS upland forest canopy is dominated by temperate deciduous early-successional species, most importantly Populus grandidentata (bigtooth aspen) and, to a lesser extent, Betula papyrifera (paper birch), with Acer saccharum (sugar maple), Acer rubrum (red maple), Fagus grandifolia (American beech), Tilia americana (basswood), Betula alleghaniensis (yellow birch), Fraxinus americana (white ash), Tsuga canadensis (eastern hemlock), Quercus rubra (northern red oak), Pinus strobus (white pine), and Pinus resinosa (red pine) existing in various fractions in the understory (and, in patches, in the canopy).
This composition is a legacy of the site's disturbance history:
The site was intensively logged and burned in the late 1800s and early 1900s, and experienced regularly recurring fires until the mid-1920s, at which point a regime of active fire suppression started that has persisted to the present day [@kilburn_1960_effects; @barnes_1981_michigan; @friedman_2005_regional_legacies_logging].
As a result, the average stand age in 2013 was 95 years [@gough_2013_sustained].
The majority of the forest that is aspen-dominated is undergoing succession to "northern hardwood" (maple, beech, basswood, birch, ash, hemlock), "upland conifer" (red and white pine), or "northern red-oak" ecotypes [@bergen_2007_observing].
Ecosystem Demography Model (version 2.2)
For this study, we used the Ecosystem Demography Model, version 2.2 (ED-2.2) to simulate forest C cycling and plant functional type (PFT) composition.
A full description of the default configuration of the model is provided by @longo_2019_ed1.
Briefly, ED-2.2 solves the energy, water, and carbon cycles separately for each of multiple "cohorts" of trees of similar composition, size, and age sharing a micro-environment and disturbance history [@moorcroft_2001_method; @medvigy_2009_mechanistic; @longo_2019_ed1].
Besides the cohort-based approximation of gap dynamics, a distinctive feature of ED-2.2 is that the equations for energy, water, and carbon fluxes are defined in terms of total energy, water, and carbon, which ensures excellent conservation of mass in long-running (multidecadal) simulations [@longo_2019_ed1].
Various versions of the ED model have been evaluated in boreal [@medvigy_2011_predicting], temperate [@medvigy_2009_mechanistic; @dietze_2014_quantitative; @raczka_2018_what], and tropical biomes [@moorcroft_2001_method; @longo_2019_ed2], and have been applied in, among others,
paleoecological studies [@rollinson_2017_emergent],
free-air CO~2~ enrichment studies [@dekauwe_2013_forest; @miller_2015_alteration],
estimation of potential carbon stocks [@hurtt_2004_beyond_potential_vegetation],
and analyses of regional vegetation-climate feedbacks [@swann_2015_future].
The official ED-2.2 source code is available at https://github.com/EDmodel/ED2, and the exact version used in this study (which includes minor revisions to accommodate the requirements of this study) can be obtained at https://github.com/ashiklom/ED2/tree/b048950.
Submodels
For our analysis of model structure, we chose to target processes of known ecological importance.
Specifically, we ran ED-2.2 with a factorial combination of the following configurations: (1) two-stream vs. multiple-scatter canopy radiative transfer models (RTMs); (2) closed canopy vs. finite crown area (a.k.a. complete vs. partial shading); and (3) static vs. plastic (varying with light level) traits.
These sub-models are summarized in Table \@ref(tab:submodels).
Both canopy RTMs in ED-2.2 resolve the full vertical radiation profile within a patch as a function of canopy structure (leaf and wood area indices, crown area, leaf angle distribution) and incident solar radiation, following the definitions in CLM 4.5 [@clm45_note].
Both models also use identical definitions for the optical properties of a single canopy layer as a function of leaf and wood optical properties, leaf angle distribution, canopy clumping, solar zenith angle, and leaf, wood, and crown area indices.
Both models represent direct (a.k.a. "beam") radiation as an exponentially decaying process, and solve for the diffuse (a.k.a. "hemispherical", "isotropic") radiation at each layer using a linear matrix equation.
The models differ in the terms of this matrix equation.
The default model (two-stream) uses a special multi-canopy solution [derived in @longo_2019_ed1] of the two-stream approximation for vegetation canopies [@dickinson_1983_land; @sellers_1985_canopy].
Meanwhile, the multiple-scatter model was derived from first principles by @zhao_2005_multiple specifically for vegetation canopies to address known biases and limitations of the two-stream approach [e.g. @wang_2003_comparison].
Compared to the default two-stream model, the multiple-scatter model slightly increases diffuse light levels in the understory.
The crown area submodel in ED-2.2 determines the nature of competition for light between cohorts.
In the default configuration (closed), the leaf area of a cohort is distributed across the entire horizontal area of a patch.
This means that taller cohorts always receive more incoming radiation---including 100% of the direct radiation---than shorter cohorts, even when the height difference is small.
Similar approaches have been shown to excessively suppress competition from sub-dominant individuals and result in unrealistically homogeneous forest patches [@fisher_2015_taking].
In the alternate configuration (finite), canopies take up only a fraction of the available horizontal area, which allows more light (including some direct radiation) to reach subdominant cohorts.
The horizontal area of crowns is determined by allometric equations from @dietze_2008_changing.
The third submodel we evaluated was trait plasticity.
In the default configuration (static), all cohorts of a given PFT have the same parameters, regardless of environmental conditions.
This ignores globally-documented intraspecific trait variability with light level and relative canopy position [@niinemets_2010_review; @keenan_2016_global].
In the alternate configuration (plastic), cohorts that receive less light increase their specific leaf area and decrease their V~cmax~, following empirical relationships from the tropics [@lloyd_2010_optimisation].
Parameter descriptions
The full list of parameters in ED-2.2 is large, with over 100 parameters per PFT.
A full sensitivity and uncertainty analysis across all of these parameters is outside the scope of this study.
Instead, we selected a subset of parameters (Table \@ref(tab:param-table)) guided by (but expanding upon) previous ED-2.2 sensitivity studies [@dietze_2014_quantitative; @raczka_2018_what; @viskari_2019_influence].
Collectively, these parameters cover a range of physiological and ecological processes relevant to C cycling and plant community composition, including
leaf physiology (light and dark reactions of photosynthesis; respiration; evapotranspiration),
carbon allocation,
plant organ stoichiometry,
canopy architecture and radiative transfer,
dispersal,
reproduction,
and mortality.
A full description of these parameters and their associated processes in ED-2.2 is provided in @longo_2019_ed1.
Plant functional types and parameterization
By default, ED-2.2 supports 17 different PFTs, which divide plant species according to growth form, photosynthetic pathway, leaf habit, biome, and successional status.
However, we limited our simulations to the four PFTs that have any appreciable presence at UMBS: Early-, mid-, and late-successional temperate deciduous ("hardwood") trees and pines.
The species comprising these PFTs are shown in Table \@ref(tab:pft-species).
For each plant functional type, we generated a distribution of parameter values via the Predictive Ecosystem Analyzer (PEcAn) trait-meta analysis [@lebauer_2013_facilitating; see also @dietze_2014_quantitative; @raczka_2018_what].
Prior distributions for this meta-analysis were largely adapted from previous ED parameter uncertainty studies [@dietze_2014_quantitative; @raczka_2018_what; @viskari_2019_influence], and are shown in detail in Supplementary Table S1.
Species trait data for this meta-analysis came from existing records in the BETY database [www.betydb.org, @lebauer_2017_betydb], as well as from publicly available records in the TRY database [www.try-db.org, @kattge_2011_try] and, specifically for leaf optical properties, from @shiklomanov_dissertation.
All parameters not described in this section were set to their ED-2.2 PFT-specific default values for all simulations.
Modeling workflow
For each factorial combination of model configurations (described above), we ran ED-2.2 once with its default parameter values, once at the posterior median estimates from the meta-analysis, and 500 times with parameters sampled randomly from the meta-analysis posterior distributions.
For the latter, the same parameter samples were used for all model structures to allow us to independently evaluate the impacts of model structure at specific parameter values.
We ran each ensemble member from January 1, 1902 through December 31, 1999 and initialized them from a "near-bare ground" condition: An equal number of seedlings of each included PFT at the minimum resolvable size.
Driving meteorological data were 6-hourly CRU-NCEP [@cruncep_pecan] combined with an annual atmospheric CO~2~ record from Law-Dome ice core [@etheridge_1998_historical] and Mauna Loa observatory [@thoning_1989_atmospheric] linearly interpolated to 3-hourly resolution.
Soil texture was set to 92% sand, 7% silt, and 1% clay, based on site-level observations [@gough_2010_wood].
The initial soil moisture profile was set to the average soil moisture profile reported in the UMBS Ameriflux ancillary data (https://ameriflux.lbl.gov/sites/siteinfo/US-UMd).
Analysis of results
First, we evaluated the impact of model structure independent of any changes to parameters by comparing predictions of average annual net primary productivity (NPP) summer (June, July, August) total leaf area index (LAI), the NPP:LAI ratio, and the effective number of PFTs across model structures while holding ED-2.2 parameters at their default values.
We calculated the effective number of PFTs ($N_{\text{PFT,eff}}$) as the inverse of the Simpson diversity index:
$$
N_{\text{PFT,eff}} = \left( \sum_i {{f_{\text{AGB,i}}}^2} \right) ^ {-1}
$$
where $f_{\text{AGB,i}}$ is the fraction of aboveground biomass for PFT $i$ at a particular time.
Second, we compared these default parameter simulations to simulations using median parameter estimates from a trait meta-analysis.
Third, we considered the full range of parameter uncertainty by looking at predictions of NPP, LAI, and composition across all parameter ensemble members.
Finally, we assessed which parameters contributed the most to parameter uncertainty for each model structure by performing a sensitivity and variance decomposition analysis following @lebauer_2013_facilitating, @dietze_2014_quantitative, and @raczka_2018_what.
Briefly, this analysis uses a spline to estimate the relationship between individual parameters and output variables and then uses these splines to predict the variance in the output variable attributable to individual parameters ["partial variance"; see @lebauer_2013_facilitating for more details].
All analyses for this work were performed using R 3.6 [@r_361].
Code and supporting data for reproducing this analysis are publicly available on the Open Science Framework (OSF) at https://osf.io/dznuf/.
Results
Impacts of structure alone
Different model structures agreed for the first half of the century-long simulation, at which point they started to diverge (Figure \@ref(fig:structure-compare-default)).
Post 1950, the combination of two-stream RTM, finite canopy radius, and plastic traits had the highest productivity, followed by the same configuration with static traits, with the remaining configurations having similar results.
The effects of changing individual model configurations were not additive---for example, enabling trait plasticity had a much stronger positive effect when paired with the two-stream canopy RTM than with the multiple-scatter canopy RTM.
Differences in productivity across model structures were related to differences in composition and vertical canopy radiation profiles over the course of the simulation (Figure S1).
Enabling the finite canopy radius increased the availability of light for understory cohorts throughout the simulation, particularly in the first half of the simulation.
Again, however, the effects of model structures were not strictly additive:
Switching from a two-stream to a multiple-scatter canopy RTM enhanced the understory light environment (albeit slightly) under a closed canopy, but the brightest understory occurred under the combination of the two-stream RTM and finite canopy.
Despite these differences, by 1980 none of the configurations was able to sustain any PFTs other than early hardwoods, and only configurations with finite canopy and two-stream RTM were able to sustain more than a single dominant cohort (two early hardwood cohorts with static traits, and three with plastic traits).
Default vs. median parameters
Parameter estimates from the trait meta-analysis were often substantially different from the ED-2.2 defaults (Figure S2).
In many cases, median parameter estimates were different entirely because of the choice of prior distributions.
However, in several important cases, parameters well-constrained by data produced values substantially different from the model defaults.
For example, default values for V~cmax~ were too low for all PFTs, and default values for SLA were too low for all PFTs except early hardwoods.
These differences in parameter values resulted in ED-2.2 runs that differed from the default runs in several important ways (Figures \@ref(fig:default-median-flux), \@ref(fig:default-median-lai)).
First, in every model configuration (except finite canopy, two-stream RTM, and static traits), gross primary productivity (GPP) reached a higher peak, and did so earlier, under median parameters (Figure \@ref(fig:default-median-flux)).
This was driven primarily by higher V~cmax~ and SLA values for most PFTs (Figure S2).
However, this increase in GPP was offset by a large increase in autotrophic respiration (RA; Figure \@ref(fig:default-median-flux)), driven primarily by higher growth and leaf respiration rates across all PFTs (notably, the ED-2.2 default is zero growth respiration for all temperate hardwood PFTs; Figure S2).
The resulting net primary productivity (NPP) had similar peak values under default and median parameters, though the temporal patterns were quite different in most model configurations (Figure \@ref(fig:default-median-flux)).
Second, the variation across model structures was much higher, with every change in model structure exerting a noticeable effect on flux variability (Figure \@ref(fig:default-median-flux)).
The most dramatic example is that, with median parameter values, simulations with the finite canopy radius produced much lower net productivity than simulations with a closed canopy in the second half of the simulation.
Again, interactions between configurations were significant; for example, under a closed canopy configuration, enabling trait plasticity caused interannual variability in NPP to decrease with a two-stream RTM and increase with a multiple-scatter RTM.
Third, even where NPP was similar between default and median runs, the underlying forest composition was different (Figure \@ref(fig:default-median-lai)).
With the default parameterization, ED-2.2 produced forests dominated by early hardwoods, with an understory of small mid- and even smaller late-hardwood trees that eventually died out.
However, with the median parameterization, ED-2.2 started out with forests dominated by mid-hardwood species that were quickly (pre-1910 to 1930, depending on configuration) replaced with pine, which remained the dominant canopy species throughout the rest of the simulation.
Exploratory analyses (results not shown) showed that the main parameters responsible for this compositional trajectory were SLA (especially the higher SLA of pine and mid-hardwood), V~cmax~, and growth respiration.
While the uncertainty in growth respiration is large, we emphasize that SLA and V~cmax~ are well constrained by data for all PFTs considered in this analysis.
Structure vs. parameter uncertainty
Net primary productivity (NPP), total leaf area index (LAI), their ratio, and PFT diversity varied with both parameter value and model structure (Figure \@ref(fig:summary-ts-plot)).
All model configurations exhibited a similar mean trajectory of rapid growth in the first 15-20 years followed by a gradual decline to a lower, stable value around 1975.
The ensemble spread within any single configuration---i.e. parameter uncertainty---was similar and large across all configurations; both NPP and LAI varied from close to zero to more than twice the mean value at any given timestep.
Within all configurations, most (~95%) ensemble members under-predicted observed present-day NPP and LAI by year 2000, though a small number of ensemble members from some configurations over-predicted one or both variables.
Peak NPP (between 1910 and 1950 for most simulations) was closer to the observed present-day NPP (though also with large ensemble spread), but LAI during this period was still too low.
The NPP to LAI ratio ("production efficiency") was higher in the closed canopy configuration than the finite canopy configuration, especially later into the simulation.
Although many simulations matched either the observed NPP or the observed LAI, only a small minority of individual runs matched both simultaneously (Figure S3).
The large parameter uncertainty in NPP and LAI also extended to species composition (Figure \@ref(fig:summary-ts-plot)).
Under all configurations, at any moment in time, most ensemble members were dominated by a single PFT, and the sustained coexistence of more than two PFTs was limited to a small set of specific ensemble members across all configurations.
However, which PFT was dominant depended strongly on parameter values, and any of the four PFTs became the dominant one in some subset of simulations (Figure S4).
Model structure and parameters interacted to determine the community assembly and dynamics of each simulation, with consequences for overall ecosystem productivity that varied by model structure (Figure \@ref(fig:lai-pft-plot)).
Parameter sets A and B both resulted in similar patterns of immediate and complete dominance by mid- and late-hardwood, respectively, followed by gradual decline to near-zero LAI.
In both of these simulations, peak productivity was higher and more sustained when the canopy was closed rather than finite, and this difference was larger under a two-stream rather than a multiple-scatter RTM.
Trait plasticity had relatively little effect on most of these two simulations, though in parameter set B, under the combination of closed canopy, multiple-scatter RTM, and plastic traits, pine overtook the dominant late hardwood PFT starting around 1975.
Both of these parameter sets also had the lowest LAI of the five selected here, with a peak LAI of 1 in A and 2 in B.
Parameter set C showed much more variety in compositional dynamics with model structure:
In all cases, mid-hardwoods began as the dominant PFT and were eventually overtaken by either pine (closed canopy and plastic traits) or early hardwood (all other cases; though with a closed canopy and static traits, the final LAI was low).
Parameter sets D and E had the highest productivity of these five sets, and were characterized by strong dominance of early hardwood and Pine, respectively.
The parameters that contributed most to predictive uncertainty in NPP during the "peak" period (1920-1950) ("partial variance") varied with model structure (Figure \@ref(fig:partial-variance-npp)).
The plant-soil water conductance (a.k.a. "water availability factor") for early hardwoods was the most or second-most important parameter for all but one model configuration (the ED-2.2 default---closed canopy, two-stream RTM, static traits---where it was not even in the top 8), and was particularly important in simulations with a finite canopy.
The growth respiration factor, particularly of Pine and, to a lesser extent, early hardwood, was also among the most influential parameters across most configurations.
Which parameters were most important, and how important they were, varied with model configuration.
For instance, the C balance ratio of mortality (mort2) of early Hardwoods was among the top four most important parameters when the canopy was closed but mattered less when the canopy was finite.
Similarly, parameters related to leaf optical properties had relatively higher importance in simulations with the two-stream model than the multiple-scatter model.
Discussion
ED-2.2 vs. reality
Our overarching objective was to investigate how model structure and parameters interact in multidecadal simulations of forest ecosystem production and community succession.
Our first question was:
Which model structural and parameter configurations reproduce observed NPP and community compositional changes?
We found that model-observation fidelity was variable, with modelled community succession generally departing from observations and modelled NPP more comparable to observed values and trends.
The community successional dynamics observed at our study site are representative of temperate forests globally, with short-lived early hardwoods initially occupying a majority of the canopy and over time gradually declining as longer-lived mid- and late-hardwood species emerge as dominant [@wolter_2002_recent; @gough_2007_legacy; @sommerfeld_2018_patterns].
At our site, Populus, a genus of early hardwoods, is abundant in the canopy during the first 50 to 75 years of succession and, upon peak maturity, declines precipitously at a rate of 40% per decade [@gough_2010_wood].
Filling this canopy void over time are mid-successional Quercus and Acer, and later successional Pinus [@scheuermann_2018_effects].
However, modelled community composition, while highly variable, overwhelmingly favored dominance by a single PFT (Figure \@ref(fig:summary-ts-plot)), with only a few parameter combinations exhibiting more realistic ecological dynamics (Figure \@ref(fig:lai-pft-plot)).
This tendency towards mono-dominance reinforces the difficulty of accurately representing plant functional diversity in vegetation models [@fisher_2015_taking].
Specific to our analysis, our results indicate that environmental filtering and competitive exclusion are potent forces in ED-2.2, and that realistic levels of coexistence and competition can only be achieved under precise combinations of the multivariate parameter space.
Unlike community succession, most model configurations produced more realistic C flux successional trajectories.
In particular, the use of median parameter values resulted in century-long C flux dynamics (Figure \@ref(fig:default-median-flux)) that mirror the trajectory of site-level observations [@gough_2007_legacy] and, more broadly, those from other temperate forests [@law_2001_carbon; @law_2003_changes; @irvine_2004_age; @tang_2009_soil; @kashian_2013_postfire].
In our study, modelled and observed C fluxes generally increased rapidly in the first quarter century before gradually declining (Figures \@ref(fig:default-median-flux) and \@ref(fig:summary-ts-plot)), a pattern that is consistent with ecological theory [@odum_1969_strategy_ecosystem_development] and global observations [@pregitzer_2004_carbon; @luyssaert_2008_old; @gough_2016_disturbance].
At odds with site-level observations, however, is the extent of simulated declines in production over time.
Rather than decreasing, long-term meteorological flux tower observations indicate net ecosystem production increased two-fold over a period of 15 years in the century-old UMBS forest, possibly because the decline of senscent early hardwoods released mid and late hardwoods poised for rapid growth [@gough_2016_disturbance; @curtis_2018_forest].
This mismatch may be driven by the aforementioned tendency of ED-2.2 towards monodominance; in the absence of an understory, there is nothing to compensate declining production among dominant trees.
When considering the coupling of composition and C fluxes, we found that the successional arc of production was similar among model configurations but improved when simulations more closely represented observed forest composition.
Most ED-2.2 runs under-predicted NPP and LAI during the peak period (1920-1950) and towards the end of the simulation (1975-1999) (Figure \@ref(fig:summary-ts-plot)).
However, of five representative runs, the two that most closely matched present-day composition at UMBS (early hardwood dominated; C and D in Figure \@ref(fig:lai-pft-plot)) also exhibited the greatest fidelity with observed NPP and LAI (Figure \@ref(fig:summary-ts-plot) and S3).
These results are encouraging, demonstrating an ecologically realistic dependency of C fluxes on canopy structure and composition [@pan_2011_large; @he_2012_relationships; @reich_2012_key].
Moreover, we show that such realistic outcomes are possible without tuning important and directly observable parameters like SLA and V~cmax~ to unrealistic values (Figure S2).
We therefore argue that, although challenging, it is not unreasonable to expect demographic vegetation models to accurately simulate both community composition and ecosystem biogeochemistry, and in fact, accuracy in the latter may depend on the former.
Drivers of uncertainty
Our second question was: What is the relative contribution of structural and parameter uncertainty to C cycling and ecological dynamics at UMBS?
We found that parameter uncertainty was by far the dominant driver of uncertainty in our simulations.
Different parameter combinations led to not only huge ranges in ensemble predictions of stand-level NPP and LAI, but also a multitude of qualitatively different ecological narratives (Figures \@ref(fig:summary-ts-plot), \@ref(fig:lai-pft-plot)).
Parameter uncertainty in demographic models therefore has broad implications for C cycle science, community ecology, and forest management.
The dominance of parameter uncertainty over structural uncertainty in our results suggests that constraining parameter uncertainty should be a higher priority than model development.
Parameter uncertainty can be reduced in two different ways:
Directly constraining parameters using observational data (for parameters that can be observed directly, such as plant functional traits) [@lebauer_2013_facilitating],
or indirectly constraining parameters through algorithmic model calibration techniques [@fer_2018_linking].
For example, leaf optical properties (i.e. reflectance, transmittance) are both measurable and highly influential [Figure \@ref(fig:partial-variance-npp); see also @viskari_2019_influence], and therefore should be a priority for new data collection.
In contrast, "growth respiration factor" is an empirical coefficient and not measurable, but highly influential [Figure \@ref(fig:partial-variance-npp); see also @dietze_2014_quantitative; @raczka_2018_what], and therefore should be constrained through calibration.
Importantly, efforts to reduce model uncertainty are most effective when they target parameters that are poorly constrained and to which the model is most sensitive.
For parameters that are already well constrained by data, the marginal value of additional measurements may be limited.
For example, ED-2.2 is highly sensitive to SLA [@dietze_2014_quantitative; @pandit_2018_optimizing], but thanks to the large number of SLA observations, we were able to estimate it very precisely (Figure S2) and, as a result, it was not a significant source of uncertainty in our analysis (Figure \@ref(fig:partial-variance-npp)).
At the same time, some poorly constrained parameters have relatively modest impacts on model outputs---for instance, rates of nonlocal dispersal and background mortality in our study (Figure S2, \@ref(fig:partial-variance-npp))---and therefore, reducing their uncertainties is not an effective way to improve model predictions.
Ultimately, this points to the value of model-data feedbacks [@dietze_2013_improving] in advancing ecological research and setting future research agendas.
Although parameter uncertainty was larger in magnitude than structural uncertainty, model structure played an important role in shaping the nature of that uncertainty.
The consequences of individual model structure changes were usually modest, but with several dramatic exceptions (Figure \@ref(fig:default-median-flux)).
Moreover, we found bi-directional interactions between structural and parameter uncertainty:
Different parameter values (especially ones far removed from ED-2.2 defaults) influenced ED-2.2 sensitivity to structural changes (Figures \@ref(fig:default-median-flux), \@ref(fig:default-median-lai)),
and changes to model structure in turn changed the relative importance of specific parameters to overall uncertainty (Figure \@ref(fig:partial-variance-npp)).
Collectively, this means that sensitivity and uncertainty analyses may not necessarily be generalizable across models, and highlights the importance of model intercomparison studies.
It also reinforces our argument that changing model structure---particularly if such changes increase model complexity---should be pursued with extreme caution.
Acknowledgments
This project was supported by NSF Award 1655095.
Cyberinfrastructure capabilities were provided by the Pacific Northwest National Laboratory (PNNL).
Additional data were provided by the University of Michigan Biological Station (UMBS) and the TRY project.
r pgbreak()
References
r pgbreak()
Tables
Table: (#tab:submodels) ED-2.2 submodel descriptions
| Name | Description | Color |
|------|-------------|-------|
| Crown model |||
| Closed (default) | Cohort crowns take up entire patch area. Competition for light based only on height. | Light |
| Finite | Cohort crown area is proportional to DBH according to PFT-specific allometry. | Dark |
| Radiative transfer model (RTM) |||
| Two-stream (default) | Two-stream approximation [@sellers_1985_canopy; @clm45_note] | Primary (red, blue) |
| Multi-scatter | Multiple-scatter approximation, following [@zhao_2005_multiple] | Secondary (green, orange) |
| Trait plasticity |||
| Static (default) | SLA and V~c,max~ are constant | Cool (blue, green) |
| Plastic | SLA increases, and V~c,max~ decreases, with light level | Warm (red, orange) |
r pgbreak()
tbl <- readr::read_csv(here(
"analysis", "data", "derived-data", "parameter-table.csv"
))
tbl <- subset(tbl, Ran == "yes", select = c(-unit_parsed, -Ran))
tbl[["ED Name"]] <- sprintf("`%s`", tbl[["ED Name"]])
cap <- paste(
"Parameter descriptions.",
"For additional information, see the",
"[ED2 model wiki](https://github.com/EDmodel/ED2/wiki/PFT-parameters)."
)
if (is_latex_output() || is_html_output()) {
kable(tbl, escape = is_latex_output(), caption = cap, booktabs = TRUE) %>%
kable_styling(latex_options = "scale_down")
} else {
# Probably DOCX output
cap <- paste("Table: (\\#tab:param-table)", cap)
flextable::flextable(tbl) %>%
flextable::autofit() %>%
flextable::fit_to_width(7.5) %>%
flextable::set_caption(cap)
}
r pgbreak()
Table: (#tab:pft-species) Plant functional type-species mappings
| Plant functional type | Species | Color |
|-----------------------|-------------------------|--------|
| Early hardwood | Betula papyrifera | Violet |
| | Populus grandidentata | |
| | Populus tremuloides | |
| Mid hardwood | Quercus rubra | Blue |
| | Acer rubrum | |
| | Acer pensylvaticum | |
| Late hardwood | Acer saccharum | Green |
| | Fagus grandifolia | |
| Pine | Pinus strobus | Yellow |
r pgbreak()
Figures
### Figure 1
structure_default_cap <- paste(
htmlfig(),
"Annual NPP and LAI from ED2 runs",
"with default parameters",
"across different configurations."
)
fig("structure-compare-default.png")
r pgbreak()
### Figure 4
default_median_flux_cap <- paste(
htmlfig(),
"ED-2.2 predictions of mean annual",
"gross (GPP) and net primary productivity (NPP)",
"and autotrophic respiration (RA) from simulations using",
"default parameter values (red) and",
"posterior median parameters from trait meta-analysis (blue)",
"across different strucures.",
"Median and default parameter values are shown in Figure S2."
)
fig("default-median-fluxes.png")
r pgbreak()
### Figure 5
default_median_lai_cap <- paste(
htmlfig(),
"ED-2.2 plot-level predictions of",
"mean growing season leaf area index (LAI)",
"by plant functional type (PFT), using",
"default parameter values (top) and",
"posterior median parameters from trait meta-analysis (bottom)",
"across different strucures."
)
fig("default-median-lai.png")
r pgbreak()
### Figure 6
summary_ts_cap <- paste(
htmlfig(),
"ED-2.2 ensemble plot-level predictions of",
"net primary productivity (NPP),",
"total leaf area index (LAI),",
"production efficiency (NPP / LAI)",
"and effective number of PFTs ($N_{\\text{PFT,eff}}$)",
"by model configuration.",
"Shading indicates the middle 50\\% (dark), middle 95\\% (medium),",
"and entire range (light) of ensemble simulations",
"(NPP/LAI results are clipped to (-5, 40) for readability,",
"as a few isolated ensemble members produced extreme values at specific",
"time points).",
"Effective number of PFTs was calculated as the inverse",
"of the Simpson diversity index with proportions expressed as",
"fraction of aboveground biomass.",
"Each colored line represents the same set of parameters across all configurations,",
"and corresponds to the letters in Figure 5.",
"Black dot with error bars is the observed mean and min/max value."
# TODO REF
)
fig("summary-ts-plot.png")
r pgbreak()
### Figure 9
lai_pft_cap <- paste(
htmlfig(),
"ED-2.2 PFT-level predictions of leaf area index (LAI) by model configuration and parameterization.",
"Each row of plots is a simulation with the same set of parameters across model configurations.",
"Letters correspond to the labels in Figure 7.",
"Note that these are predictions of total leaf area index summed across all cohorts of the same PFT,",
"meaning that higher LAI in one PFT does not necessarily indicate dominance, size, etc."
)
fig("lai-pft.png")
r pgbreak()
### Figure 10
partial_variance_npp_cap <- paste(
htmlfig(),
"Contributions of individual parameters to predicted uncertainty in",
"ED-2.2 predictions of average net primary productivity (NPP)",
"over the period 1920-1950",
"by model type.",
"Each panel shows the top 8 parameter values by their fractional contribution",
'to parameter uncertainty ("partial variance").',
"Input parameter distributions are shown in Figure 3.",
"Note that x-axis scales vary across panels."
)
fig("partial-variance-npp.png")
r pgbreak()
ashiklom/fortebaseline documentation built on May 9, 2020, 1:56 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
library(knitr) library(here) library(magrittr) library(kableExtra) fig <- function(...) { knitr::include_graphics(here::here("analysis", "figures", ...)) } out_format <- knitr::opts_knit$get("rmarkdown.pandoc.to") pgbreak <- function() { if (grepl("docx", out_format)) { return("##### pagebreak") } else if (is_latex_output()) { return("\\newpage") } } opts_chunk$set( collapse = TRUE, warning = FALSE, message = FALSE, echo = FALSE, comment = "#>", out.width = "\\linewidth" ## fig.path = "../figures/" ) ## library(fortebaseline) # Or use devtools::load_all('.', quiet = T) if your code is in script files, rather than as functions in the `/R` diretory n <- 0 mdcap <- function(text) { n <<- n + 1 out <- sprintf("**Figure %d**: %s", n, text) if (grepl("gfm", out_format, ignore.case = TRUE)) { out } else { invisible(NULL) } } m <- 0 htmlfig <- function() { m <<- m + 1 out <- paste0("**Figure ", m, "**:") if (knitr::is_html_output()) { out } else { invisible(NULL) } }
Introduction
Human land use patterns over the last 500 years, and particularly over the last century, have led to decline in primary forests across most of the Earth's surface [@hurtt_2011_harmonization]. However, in many places, subsequent changes in land use (e.g. abandonment of agriculture and pasture) have caused much of this area to regrow as secondary forest [@hurtt_2011_harmonization]. Understanding the drivers behind this regrowth, and the ability to predict its future course, are therefore essential to shaping global policy around climate, biodiversity, and natural resource management.
The temperate forests of the Upper Midwest are a model system for studying secondary forest ecology and biogeochemistry. These forests have been a net carbon sink for most of the past century, as fast-growing early successional species established and thrived following near-complete deforestation around the end of the 19th Century [@birdsey_2006_forest; @williams_2012_carbon]. More recently, many Upper Midwest forests are undergoing a shift in community composition from early-successional to mid-successional species, and the ecological and biogeochemical consequences of this transition are not well-established [@gough_2016_disturbance]. In addition, these forests face a multitude of direct and indirect pressures from both biotic and abiotic sources [@shifley_2016_future], including longer and more severe droughts [@gustafson_2016_implications; @swanston_2017_vulnerability; @lienard_2016_us] and non-native insects and pathogens [@lovett_2016_nonnative]. Collectively, these processes have culminated in forests that are different from their pre-settlement climax counterparts [@stearns_2002_one; @wolter_2002_recent; @thompson_2013_four], and our ability to make predictions about such non-analog environments based on the past is limited.
More reliable predictions can in theory be obtained by using dynamic vegetation models that explicitly represent processes involved in forest growth and mortality. Vegetation models fall along a spectrum of complexity, with trade-offs between realism and process fidelity on one hand and computational demand and data requirements on the other. On one side are "big leaf" models---such as PNET [@aber_1992_pnet], SiBCASA [@schaefer_2008_sibcasa], and Biome-BGC [@running_1993_generalization]---in which the vegetation at a particular location is represented by a single average "plant". On the other are "individual-based" (or "gap") models---such as LANDIS [@mladenoff_2004_landis; @scheller_2007_design] and UVAFME [@shuman_2015_forest]---that explicitly simulate competition between multiple individuals [@shugart_2015_computer]. In the middle are models that approximate the emergent outcomes of competitive interactions among different individuals without explicitly simulating each individual [@purves_2008_predictive]. One example of this approach is the Ecosystem Demography model [ED; @moorcroft_2001_method; @medvigy_2011_predicting; @longo_2019_ed1], which is the focus of this study.
Regardless of vegetation model complexity, their projections are inherently uncertain. This uncertainty comes from many different sources, of which we consider two in this study: Structural uncertainty refers to which process are included and how they are represented. Model intercomparison studies have shown the importance of different model structures to projections of overall land carbon sequestration [@friedlingstein_2006_climate; @friedlingstein_2014_uncertainties; @lovenduski_2017_reducing], response to CO~2~ fertilization [@zaehle_2014_evaluation; @medlyn_2015_using; @walker_2015_predicting], and soil C sequestration [@sulman_2018_multiple]. This structural uncertainty has been attributed to differences in model representations of key processes, including canopy radiative transfer [@fisher_2017_vegetation], soil biogeochemistry [@wieder_2017_carbon; @sulman_2018_multiple], stomatal conductance [@franks_2018_comparing], and photosynthesis [@rogers_2016_roadmap]. However, the extent to which specific processes contribute to model uncertainty is difficult to evaluate from model intercomparisons because models differ from each other in too many ways, meaning that structural effects are confounded with other uncertainty sources. At the same time, parameter uncertainty arises because of natural variability and imperfect calibration in process representations. Previous analyses of the ED model in particular have highlighted the importance of parameter uncertainty (and reduction thereof) to predictions of primary productivity and composition [@dietze_2014_quantitative; @raczka_2018_what; @viskari_2019_influence]. While many modeling studies have looked at structural and parameter uncertainty independently, comparative analysis of contributions of both uncertainty sources are rare in the ecosystem modeling literature.
Light availability and distribution drives demographic variation in plant C fixation and community structural changes through succession. Observations show that competition for light drives niche differentiation [@silvertown_2004_plant] and consequently successional change [@canham_1990_light; @bazzaz_1993_successional], particularly in closed canopies. Prior ED versions were limited in their capacity to sustain multiple plant functional types during stand development due in part to unrealistically aggressive shading of the understory and a lack of consideration for tree crown area, such that the marginal height advantage of a vegetation cohort resulted in a disproportionally large growth advantage [@fisher_2010_assessing]. Alternative approaches that accommodate horizontal competition for light [e.g. perfect plasticity assumption, @strigul_2008_scaling; @weng_2015_scaling] reduce the relative impact of height on light availability and growth, and may therefore promote competition [@fisher_2017_vegetation]. In addition, trait plasticity may permit individuals to acclimate and species to adapt to changing light environments in ways that improve light use efficiency [@niinemets_1998_analysis; @lloyd_2010_optimisation; @keenan_2016_global]. To this end, we focused our structural uncertainty analysis on three processes related to competition for- and utilization efficiency of light: Representation of horizontal competition; representation of canopy radiative transfer (especially of diffuse radiation); and trait plasticity with light level.
This study focuses on the interaction between parameter and structural uncertainty. Our study is organized around the following guiding questions: (1) Which processes related to light utilization are most important to accurately model community succession and C cycling in Upper Midwest forests? (2) What are the relative contributions of parameter uncertainty (data limitation) vs. structural uncertainty (theoretical limitation)? We ran ensemble simulations of the Ecosystem Demography (ED-2.2) model with a factorial combination of submodels related to canopy radiative transfer formulation, horizontal competition, and trait plasticity. We then used sensitivity and variance decomposition analyses to evaluate the contribution of parameter uncertainty for each model configuration, exploring the ecological implications of the results in the context of ongoing global change.
Methods
Site description
We performed this study at the University of Michigan Biological Station (UMBS; Ameriflux site US-UMd, 45.5625$^{\circ}$, -84.6975$^{\circ}$), located in Northern Lower Michigan, USA. The area surrounding the research station is 87% well-drained upland forest and 13% wetland [@bergen_2007_observing]; the focus of this study is on the former. The landscape geography of the UMBS upland forest is 20.4% moraine, 37.8% high outwash plain, 31.3% low outwash plain, 5.7% lake-margin terrace, 3.6% ice-contact, and 1.2% lowland glacial lake [@bergen_2007_observing]. The UMBS upland forest canopy is dominated by temperate deciduous early-successional species, most importantly Populus grandidentata (bigtooth aspen) and, to a lesser extent, Betula papyrifera (paper birch), with Acer saccharum (sugar maple), Acer rubrum (red maple), Fagus grandifolia (American beech), Tilia americana (basswood), Betula alleghaniensis (yellow birch), Fraxinus americana (white ash), Tsuga canadensis (eastern hemlock), Quercus rubra (northern red oak), Pinus strobus (white pine), and Pinus resinosa (red pine) existing in various fractions in the understory (and, in patches, in the canopy). This composition is a legacy of the site's disturbance history: The site was intensively logged and burned in the late 1800s and early 1900s, and experienced regularly recurring fires until the mid-1920s, at which point a regime of active fire suppression started that has persisted to the present day [@kilburn_1960_effects; @barnes_1981_michigan; @friedman_2005_regional_legacies_logging]. As a result, the average stand age in 2013 was 95 years [@gough_2013_sustained]. The majority of the forest that is aspen-dominated is undergoing succession to "northern hardwood" (maple, beech, basswood, birch, ash, hemlock), "upland conifer" (red and white pine), or "northern red-oak" ecotypes [@bergen_2007_observing].
Ecosystem Demography Model (version 2.2)
For this study, we used the Ecosystem Demography Model, version 2.2 (ED-2.2) to simulate forest C cycling and plant functional type (PFT) composition. A full description of the default configuration of the model is provided by @longo_2019_ed1. Briefly, ED-2.2 solves the energy, water, and carbon cycles separately for each of multiple "cohorts" of trees of similar composition, size, and age sharing a micro-environment and disturbance history [@moorcroft_2001_method; @medvigy_2009_mechanistic; @longo_2019_ed1]. Besides the cohort-based approximation of gap dynamics, a distinctive feature of ED-2.2 is that the equations for energy, water, and carbon fluxes are defined in terms of total energy, water, and carbon, which ensures excellent conservation of mass in long-running (multidecadal) simulations [@longo_2019_ed1]. Various versions of the ED model have been evaluated in boreal [@medvigy_2011_predicting], temperate [@medvigy_2009_mechanistic; @dietze_2014_quantitative; @raczka_2018_what], and tropical biomes [@moorcroft_2001_method; @longo_2019_ed2], and have been applied in, among others, paleoecological studies [@rollinson_2017_emergent], free-air CO~2~ enrichment studies [@dekauwe_2013_forest; @miller_2015_alteration], estimation of potential carbon stocks [@hurtt_2004_beyond_potential_vegetation], and analyses of regional vegetation-climate feedbacks [@swann_2015_future].
The official ED-2.2 source code is available at https://github.com/EDmodel/ED2, and the exact version used in this study (which includes minor revisions to accommodate the requirements of this study) can be obtained at https://github.com/ashiklom/ED2/tree/b048950.
Submodels
For our analysis of model structure, we chose to target processes of known ecological importance. Specifically, we ran ED-2.2 with a factorial combination of the following configurations: (1) two-stream vs. multiple-scatter canopy radiative transfer models (RTMs); (2) closed canopy vs. finite crown area (a.k.a. complete vs. partial shading); and (3) static vs. plastic (varying with light level) traits. These sub-models are summarized in Table \@ref(tab:submodels).
Both canopy RTMs in ED-2.2 resolve the full vertical radiation profile within a patch as a function of canopy structure (leaf and wood area indices, crown area, leaf angle distribution) and incident solar radiation, following the definitions in CLM 4.5 [@clm45_note]. Both models also use identical definitions for the optical properties of a single canopy layer as a function of leaf and wood optical properties, leaf angle distribution, canopy clumping, solar zenith angle, and leaf, wood, and crown area indices. Both models represent direct (a.k.a. "beam") radiation as an exponentially decaying process, and solve for the diffuse (a.k.a. "hemispherical", "isotropic") radiation at each layer using a linear matrix equation. The models differ in the terms of this matrix equation. The default model (two-stream) uses a special multi-canopy solution [derived in @longo_2019_ed1] of the two-stream approximation for vegetation canopies [@dickinson_1983_land; @sellers_1985_canopy]. Meanwhile, the multiple-scatter model was derived from first principles by @zhao_2005_multiple specifically for vegetation canopies to address known biases and limitations of the two-stream approach [e.g. @wang_2003_comparison]. Compared to the default two-stream model, the multiple-scatter model slightly increases diffuse light levels in the understory.
The crown area submodel in ED-2.2 determines the nature of competition for light between cohorts. In the default configuration (closed), the leaf area of a cohort is distributed across the entire horizontal area of a patch. This means that taller cohorts always receive more incoming radiation---including 100% of the direct radiation---than shorter cohorts, even when the height difference is small. Similar approaches have been shown to excessively suppress competition from sub-dominant individuals and result in unrealistically homogeneous forest patches [@fisher_2015_taking]. In the alternate configuration (finite), canopies take up only a fraction of the available horizontal area, which allows more light (including some direct radiation) to reach subdominant cohorts. The horizontal area of crowns is determined by allometric equations from @dietze_2008_changing.
The third submodel we evaluated was trait plasticity. In the default configuration (static), all cohorts of a given PFT have the same parameters, regardless of environmental conditions. This ignores globally-documented intraspecific trait variability with light level and relative canopy position [@niinemets_2010_review; @keenan_2016_global]. In the alternate configuration (plastic), cohorts that receive less light increase their specific leaf area and decrease their V~cmax~, following empirical relationships from the tropics [@lloyd_2010_optimisation].
Parameter descriptions
The full list of parameters in ED-2.2 is large, with over 100 parameters per PFT. A full sensitivity and uncertainty analysis across all of these parameters is outside the scope of this study. Instead, we selected a subset of parameters (Table \@ref(tab:param-table)) guided by (but expanding upon) previous ED-2.2 sensitivity studies [@dietze_2014_quantitative; @raczka_2018_what; @viskari_2019_influence]. Collectively, these parameters cover a range of physiological and ecological processes relevant to C cycling and plant community composition, including leaf physiology (light and dark reactions of photosynthesis; respiration; evapotranspiration), carbon allocation, plant organ stoichiometry, canopy architecture and radiative transfer, dispersal, reproduction, and mortality. A full description of these parameters and their associated processes in ED-2.2 is provided in @longo_2019_ed1.
Plant functional types and parameterization
By default, ED-2.2 supports 17 different PFTs, which divide plant species according to growth form, photosynthetic pathway, leaf habit, biome, and successional status. However, we limited our simulations to the four PFTs that have any appreciable presence at UMBS: Early-, mid-, and late-successional temperate deciduous ("hardwood") trees and pines. The species comprising these PFTs are shown in Table \@ref(tab:pft-species).
For each plant functional type, we generated a distribution of parameter values via the Predictive Ecosystem Analyzer (PEcAn) trait-meta analysis [@lebauer_2013_facilitating; see also @dietze_2014_quantitative; @raczka_2018_what]. Prior distributions for this meta-analysis were largely adapted from previous ED parameter uncertainty studies [@dietze_2014_quantitative; @raczka_2018_what; @viskari_2019_influence], and are shown in detail in Supplementary Table S1. Species trait data for this meta-analysis came from existing records in the BETY database [www.betydb.org, @lebauer_2017_betydb], as well as from publicly available records in the TRY database [www.try-db.org, @kattge_2011_try] and, specifically for leaf optical properties, from @shiklomanov_dissertation. All parameters not described in this section were set to their ED-2.2 PFT-specific default values for all simulations.
Modeling workflow
For each factorial combination of model configurations (described above), we ran ED-2.2 once with its default parameter values, once at the posterior median estimates from the meta-analysis, and 500 times with parameters sampled randomly from the meta-analysis posterior distributions. For the latter, the same parameter samples were used for all model structures to allow us to independently evaluate the impacts of model structure at specific parameter values. We ran each ensemble member from January 1, 1902 through December 31, 1999 and initialized them from a "near-bare ground" condition: An equal number of seedlings of each included PFT at the minimum resolvable size. Driving meteorological data were 6-hourly CRU-NCEP [@cruncep_pecan] combined with an annual atmospheric CO~2~ record from Law-Dome ice core [@etheridge_1998_historical] and Mauna Loa observatory [@thoning_1989_atmospheric] linearly interpolated to 3-hourly resolution. Soil texture was set to 92% sand, 7% silt, and 1% clay, based on site-level observations [@gough_2010_wood]. The initial soil moisture profile was set to the average soil moisture profile reported in the UMBS Ameriflux ancillary data (https://ameriflux.lbl.gov/sites/siteinfo/US-UMd).
Analysis of results
First, we evaluated the impact of model structure independent of any changes to parameters by comparing predictions of average annual net primary productivity (NPP) summer (June, July, August) total leaf area index (LAI), the NPP:LAI ratio, and the effective number of PFTs across model structures while holding ED-2.2 parameters at their default values. We calculated the effective number of PFTs ($N_{\text{PFT,eff}}$) as the inverse of the Simpson diversity index:
$$ N_{\text{PFT,eff}} = \left( \sum_i {{f_{\text{AGB,i}}}^2} \right) ^ {-1} $$
where $f_{\text{AGB,i}}$ is the fraction of aboveground biomass for PFT $i$ at a particular time.
Second, we compared these default parameter simulations to simulations using median parameter estimates from a trait meta-analysis. Third, we considered the full range of parameter uncertainty by looking at predictions of NPP, LAI, and composition across all parameter ensemble members. Finally, we assessed which parameters contributed the most to parameter uncertainty for each model structure by performing a sensitivity and variance decomposition analysis following @lebauer_2013_facilitating, @dietze_2014_quantitative, and @raczka_2018_what. Briefly, this analysis uses a spline to estimate the relationship between individual parameters and output variables and then uses these splines to predict the variance in the output variable attributable to individual parameters ["partial variance"; see @lebauer_2013_facilitating for more details].
All analyses for this work were performed using R 3.6 [@r_361]. Code and supporting data for reproducing this analysis are publicly available on the Open Science Framework (OSF) at https://osf.io/dznuf/.
Results
Impacts of structure alone
Different model structures agreed for the first half of the century-long simulation, at which point they started to diverge (Figure \@ref(fig:structure-compare-default)). Post 1950, the combination of two-stream RTM, finite canopy radius, and plastic traits had the highest productivity, followed by the same configuration with static traits, with the remaining configurations having similar results. The effects of changing individual model configurations were not additive---for example, enabling trait plasticity had a much stronger positive effect when paired with the two-stream canopy RTM than with the multiple-scatter canopy RTM.
Differences in productivity across model structures were related to differences in composition and vertical canopy radiation profiles over the course of the simulation (Figure S1). Enabling the finite canopy radius increased the availability of light for understory cohorts throughout the simulation, particularly in the first half of the simulation. Again, however, the effects of model structures were not strictly additive: Switching from a two-stream to a multiple-scatter canopy RTM enhanced the understory light environment (albeit slightly) under a closed canopy, but the brightest understory occurred under the combination of the two-stream RTM and finite canopy. Despite these differences, by 1980 none of the configurations was able to sustain any PFTs other than early hardwoods, and only configurations with finite canopy and two-stream RTM were able to sustain more than a single dominant cohort (two early hardwood cohorts with static traits, and three with plastic traits).
Default vs. median parameters
Parameter estimates from the trait meta-analysis were often substantially different from the ED-2.2 defaults (Figure S2). In many cases, median parameter estimates were different entirely because of the choice of prior distributions. However, in several important cases, parameters well-constrained by data produced values substantially different from the model defaults. For example, default values for V~cmax~ were too low for all PFTs, and default values for SLA were too low for all PFTs except early hardwoods.
These differences in parameter values resulted in ED-2.2 runs that differed from the default runs in several important ways (Figures \@ref(fig:default-median-flux), \@ref(fig:default-median-lai)). First, in every model configuration (except finite canopy, two-stream RTM, and static traits), gross primary productivity (GPP) reached a higher peak, and did so earlier, under median parameters (Figure \@ref(fig:default-median-flux)). This was driven primarily by higher V~cmax~ and SLA values for most PFTs (Figure S2). However, this increase in GPP was offset by a large increase in autotrophic respiration (RA; Figure \@ref(fig:default-median-flux)), driven primarily by higher growth and leaf respiration rates across all PFTs (notably, the ED-2.2 default is zero growth respiration for all temperate hardwood PFTs; Figure S2). The resulting net primary productivity (NPP) had similar peak values under default and median parameters, though the temporal patterns were quite different in most model configurations (Figure \@ref(fig:default-median-flux)).
Second, the variation across model structures was much higher, with every change in model structure exerting a noticeable effect on flux variability (Figure \@ref(fig:default-median-flux)). The most dramatic example is that, with median parameter values, simulations with the finite canopy radius produced much lower net productivity than simulations with a closed canopy in the second half of the simulation. Again, interactions between configurations were significant; for example, under a closed canopy configuration, enabling trait plasticity caused interannual variability in NPP to decrease with a two-stream RTM and increase with a multiple-scatter RTM.
Third, even where NPP was similar between default and median runs, the underlying forest composition was different (Figure \@ref(fig:default-median-lai)). With the default parameterization, ED-2.2 produced forests dominated by early hardwoods, with an understory of small mid- and even smaller late-hardwood trees that eventually died out. However, with the median parameterization, ED-2.2 started out with forests dominated by mid-hardwood species that were quickly (pre-1910 to 1930, depending on configuration) replaced with pine, which remained the dominant canopy species throughout the rest of the simulation. Exploratory analyses (results not shown) showed that the main parameters responsible for this compositional trajectory were SLA (especially the higher SLA of pine and mid-hardwood), V~cmax~, and growth respiration. While the uncertainty in growth respiration is large, we emphasize that SLA and V~cmax~ are well constrained by data for all PFTs considered in this analysis.
Structure vs. parameter uncertainty
Net primary productivity (NPP), total leaf area index (LAI), their ratio, and PFT diversity varied with both parameter value and model structure (Figure \@ref(fig:summary-ts-plot)). All model configurations exhibited a similar mean trajectory of rapid growth in the first 15-20 years followed by a gradual decline to a lower, stable value around 1975. The ensemble spread within any single configuration---i.e. parameter uncertainty---was similar and large across all configurations; both NPP and LAI varied from close to zero to more than twice the mean value at any given timestep. Within all configurations, most (~95%) ensemble members under-predicted observed present-day NPP and LAI by year 2000, though a small number of ensemble members from some configurations over-predicted one or both variables. Peak NPP (between 1910 and 1950 for most simulations) was closer to the observed present-day NPP (though also with large ensemble spread), but LAI during this period was still too low. The NPP to LAI ratio ("production efficiency") was higher in the closed canopy configuration than the finite canopy configuration, especially later into the simulation. Although many simulations matched either the observed NPP or the observed LAI, only a small minority of individual runs matched both simultaneously (Figure S3).
The large parameter uncertainty in NPP and LAI also extended to species composition (Figure \@ref(fig:summary-ts-plot)). Under all configurations, at any moment in time, most ensemble members were dominated by a single PFT, and the sustained coexistence of more than two PFTs was limited to a small set of specific ensemble members across all configurations. However, which PFT was dominant depended strongly on parameter values, and any of the four PFTs became the dominant one in some subset of simulations (Figure S4).
Model structure and parameters interacted to determine the community assembly and dynamics of each simulation, with consequences for overall ecosystem productivity that varied by model structure (Figure \@ref(fig:lai-pft-plot)). Parameter sets A and B both resulted in similar patterns of immediate and complete dominance by mid- and late-hardwood, respectively, followed by gradual decline to near-zero LAI. In both of these simulations, peak productivity was higher and more sustained when the canopy was closed rather than finite, and this difference was larger under a two-stream rather than a multiple-scatter RTM. Trait plasticity had relatively little effect on most of these two simulations, though in parameter set B, under the combination of closed canopy, multiple-scatter RTM, and plastic traits, pine overtook the dominant late hardwood PFT starting around 1975. Both of these parameter sets also had the lowest LAI of the five selected here, with a peak LAI of 1 in A and 2 in B. Parameter set C showed much more variety in compositional dynamics with model structure: In all cases, mid-hardwoods began as the dominant PFT and were eventually overtaken by either pine (closed canopy and plastic traits) or early hardwood (all other cases; though with a closed canopy and static traits, the final LAI was low). Parameter sets D and E had the highest productivity of these five sets, and were characterized by strong dominance of early hardwood and Pine, respectively.
The parameters that contributed most to predictive uncertainty in NPP during the "peak" period (1920-1950) ("partial variance") varied with model structure (Figure \@ref(fig:partial-variance-npp)).
The plant-soil water conductance (a.k.a. "water availability factor") for early hardwoods was the most or second-most important parameter for all but one model configuration (the ED-2.2 default---closed canopy, two-stream RTM, static traits---where it was not even in the top 8), and was particularly important in simulations with a finite canopy.
The growth respiration factor, particularly of Pine and, to a lesser extent, early hardwood, was also among the most influential parameters across most configurations.
Which parameters were most important, and how important they were, varied with model configuration.
For instance, the C balance ratio of mortality (mort2) of early Hardwoods was among the top four most important parameters when the canopy was closed but mattered less when the canopy was finite.
Similarly, parameters related to leaf optical properties had relatively higher importance in simulations with the two-stream model than the multiple-scatter model.
Discussion
ED-2.2 vs. reality
Our overarching objective was to investigate how model structure and parameters interact in multidecadal simulations of forest ecosystem production and community succession. Our first question was: Which model structural and parameter configurations reproduce observed NPP and community compositional changes? We found that model-observation fidelity was variable, with modelled community succession generally departing from observations and modelled NPP more comparable to observed values and trends.
The community successional dynamics observed at our study site are representative of temperate forests globally, with short-lived early hardwoods initially occupying a majority of the canopy and over time gradually declining as longer-lived mid- and late-hardwood species emerge as dominant [@wolter_2002_recent; @gough_2007_legacy; @sommerfeld_2018_patterns]. At our site, Populus, a genus of early hardwoods, is abundant in the canopy during the first 50 to 75 years of succession and, upon peak maturity, declines precipitously at a rate of 40% per decade [@gough_2010_wood]. Filling this canopy void over time are mid-successional Quercus and Acer, and later successional Pinus [@scheuermann_2018_effects]. However, modelled community composition, while highly variable, overwhelmingly favored dominance by a single PFT (Figure \@ref(fig:summary-ts-plot)), with only a few parameter combinations exhibiting more realistic ecological dynamics (Figure \@ref(fig:lai-pft-plot)). This tendency towards mono-dominance reinforces the difficulty of accurately representing plant functional diversity in vegetation models [@fisher_2015_taking]. Specific to our analysis, our results indicate that environmental filtering and competitive exclusion are potent forces in ED-2.2, and that realistic levels of coexistence and competition can only be achieved under precise combinations of the multivariate parameter space.
Unlike community succession, most model configurations produced more realistic C flux successional trajectories. In particular, the use of median parameter values resulted in century-long C flux dynamics (Figure \@ref(fig:default-median-flux)) that mirror the trajectory of site-level observations [@gough_2007_legacy] and, more broadly, those from other temperate forests [@law_2001_carbon; @law_2003_changes; @irvine_2004_age; @tang_2009_soil; @kashian_2013_postfire]. In our study, modelled and observed C fluxes generally increased rapidly in the first quarter century before gradually declining (Figures \@ref(fig:default-median-flux) and \@ref(fig:summary-ts-plot)), a pattern that is consistent with ecological theory [@odum_1969_strategy_ecosystem_development] and global observations [@pregitzer_2004_carbon; @luyssaert_2008_old; @gough_2016_disturbance]. At odds with site-level observations, however, is the extent of simulated declines in production over time. Rather than decreasing, long-term meteorological flux tower observations indicate net ecosystem production increased two-fold over a period of 15 years in the century-old UMBS forest, possibly because the decline of senscent early hardwoods released mid and late hardwoods poised for rapid growth [@gough_2016_disturbance; @curtis_2018_forest]. This mismatch may be driven by the aforementioned tendency of ED-2.2 towards monodominance; in the absence of an understory, there is nothing to compensate declining production among dominant trees.
When considering the coupling of composition and C fluxes, we found that the successional arc of production was similar among model configurations but improved when simulations more closely represented observed forest composition. Most ED-2.2 runs under-predicted NPP and LAI during the peak period (1920-1950) and towards the end of the simulation (1975-1999) (Figure \@ref(fig:summary-ts-plot)). However, of five representative runs, the two that most closely matched present-day composition at UMBS (early hardwood dominated; C and D in Figure \@ref(fig:lai-pft-plot)) also exhibited the greatest fidelity with observed NPP and LAI (Figure \@ref(fig:summary-ts-plot) and S3). These results are encouraging, demonstrating an ecologically realistic dependency of C fluxes on canopy structure and composition [@pan_2011_large; @he_2012_relationships; @reich_2012_key]. Moreover, we show that such realistic outcomes are possible without tuning important and directly observable parameters like SLA and V~cmax~ to unrealistic values (Figure S2). We therefore argue that, although challenging, it is not unreasonable to expect demographic vegetation models to accurately simulate both community composition and ecosystem biogeochemistry, and in fact, accuracy in the latter may depend on the former.
Drivers of uncertainty
Our second question was: What is the relative contribution of structural and parameter uncertainty to C cycling and ecological dynamics at UMBS? We found that parameter uncertainty was by far the dominant driver of uncertainty in our simulations. Different parameter combinations led to not only huge ranges in ensemble predictions of stand-level NPP and LAI, but also a multitude of qualitatively different ecological narratives (Figures \@ref(fig:summary-ts-plot), \@ref(fig:lai-pft-plot)). Parameter uncertainty in demographic models therefore has broad implications for C cycle science, community ecology, and forest management.
The dominance of parameter uncertainty over structural uncertainty in our results suggests that constraining parameter uncertainty should be a higher priority than model development. Parameter uncertainty can be reduced in two different ways: Directly constraining parameters using observational data (for parameters that can be observed directly, such as plant functional traits) [@lebauer_2013_facilitating], or indirectly constraining parameters through algorithmic model calibration techniques [@fer_2018_linking]. For example, leaf optical properties (i.e. reflectance, transmittance) are both measurable and highly influential [Figure \@ref(fig:partial-variance-npp); see also @viskari_2019_influence], and therefore should be a priority for new data collection. In contrast, "growth respiration factor" is an empirical coefficient and not measurable, but highly influential [Figure \@ref(fig:partial-variance-npp); see also @dietze_2014_quantitative; @raczka_2018_what], and therefore should be constrained through calibration.
Importantly, efforts to reduce model uncertainty are most effective when they target parameters that are poorly constrained and to which the model is most sensitive. For parameters that are already well constrained by data, the marginal value of additional measurements may be limited. For example, ED-2.2 is highly sensitive to SLA [@dietze_2014_quantitative; @pandit_2018_optimizing], but thanks to the large number of SLA observations, we were able to estimate it very precisely (Figure S2) and, as a result, it was not a significant source of uncertainty in our analysis (Figure \@ref(fig:partial-variance-npp)). At the same time, some poorly constrained parameters have relatively modest impacts on model outputs---for instance, rates of nonlocal dispersal and background mortality in our study (Figure S2, \@ref(fig:partial-variance-npp))---and therefore, reducing their uncertainties is not an effective way to improve model predictions. Ultimately, this points to the value of model-data feedbacks [@dietze_2013_improving] in advancing ecological research and setting future research agendas.
Although parameter uncertainty was larger in magnitude than structural uncertainty, model structure played an important role in shaping the nature of that uncertainty. The consequences of individual model structure changes were usually modest, but with several dramatic exceptions (Figure \@ref(fig:default-median-flux)). Moreover, we found bi-directional interactions between structural and parameter uncertainty: Different parameter values (especially ones far removed from ED-2.2 defaults) influenced ED-2.2 sensitivity to structural changes (Figures \@ref(fig:default-median-flux), \@ref(fig:default-median-lai)), and changes to model structure in turn changed the relative importance of specific parameters to overall uncertainty (Figure \@ref(fig:partial-variance-npp)). Collectively, this means that sensitivity and uncertainty analyses may not necessarily be generalizable across models, and highlights the importance of model intercomparison studies. It also reinforces our argument that changing model structure---particularly if such changes increase model complexity---should be pursued with extreme caution.
Acknowledgments
This project was supported by NSF Award 1655095. Cyberinfrastructure capabilities were provided by the Pacific Northwest National Laboratory (PNNL). Additional data were provided by the University of Michigan Biological Station (UMBS) and the TRY project.
r pgbreak()
References
r pgbreak()
Tables
Table: (#tab:submodels) ED-2.2 submodel descriptions
| Name | Description | Color | |------|-------------|-------| | Crown model ||| | Closed (default) | Cohort crowns take up entire patch area. Competition for light based only on height. | Light | | Finite | Cohort crown area is proportional to DBH according to PFT-specific allometry. | Dark | | Radiative transfer model (RTM) ||| | Two-stream (default) | Two-stream approximation [@sellers_1985_canopy; @clm45_note] | Primary (red, blue) | | Multi-scatter | Multiple-scatter approximation, following [@zhao_2005_multiple] | Secondary (green, orange) | | Trait plasticity ||| | Static (default) | SLA and V~c,max~ are constant | Cool (blue, green) | | Plastic | SLA increases, and V~c,max~ decreases, with light level | Warm (red, orange) |
r pgbreak()
tbl <- readr::read_csv(here( "analysis", "data", "derived-data", "parameter-table.csv" )) tbl <- subset(tbl, Ran == "yes", select = c(-unit_parsed, -Ran)) tbl[["ED Name"]] <- sprintf("`%s`", tbl[["ED Name"]]) cap <- paste( "Parameter descriptions.", "For additional information, see the", "[ED2 model wiki](https://github.com/EDmodel/ED2/wiki/PFT-parameters)." ) if (is_latex_output() || is_html_output()) { kable(tbl, escape = is_latex_output(), caption = cap, booktabs = TRUE) %>% kable_styling(latex_options = "scale_down") } else { # Probably DOCX output cap <- paste("Table: (\\#tab:param-table)", cap) flextable::flextable(tbl) %>% flextable::autofit() %>% flextable::fit_to_width(7.5) %>% flextable::set_caption(cap) }
r pgbreak()
Table: (#tab:pft-species) Plant functional type-species mappings
| Plant functional type | Species | Color | |-----------------------|-------------------------|--------| | Early hardwood | Betula papyrifera | Violet | | | Populus grandidentata | | | | Populus tremuloides | | | Mid hardwood | Quercus rubra | Blue | | | Acer rubrum | | | | Acer pensylvaticum | | | Late hardwood | Acer saccharum | Green | | | Fagus grandifolia | | | Pine | Pinus strobus | Yellow |
r pgbreak()
Figures
### Figure 1 structure_default_cap <- paste( htmlfig(), "Annual NPP and LAI from ED2 runs", "with default parameters", "across different configurations." ) fig("structure-compare-default.png")
r pgbreak()
### Figure 4 default_median_flux_cap <- paste( htmlfig(), "ED-2.2 predictions of mean annual", "gross (GPP) and net primary productivity (NPP)", "and autotrophic respiration (RA) from simulations using", "default parameter values (red) and", "posterior median parameters from trait meta-analysis (blue)", "across different strucures.", "Median and default parameter values are shown in Figure S2." ) fig("default-median-fluxes.png")
r pgbreak()
### Figure 5 default_median_lai_cap <- paste( htmlfig(), "ED-2.2 plot-level predictions of", "mean growing season leaf area index (LAI)", "by plant functional type (PFT), using", "default parameter values (top) and", "posterior median parameters from trait meta-analysis (bottom)", "across different strucures." ) fig("default-median-lai.png")
r pgbreak()
### Figure 6 summary_ts_cap <- paste( htmlfig(), "ED-2.2 ensemble plot-level predictions of", "net primary productivity (NPP),", "total leaf area index (LAI),", "production efficiency (NPP / LAI)", "and effective number of PFTs ($N_{\\text{PFT,eff}}$)", "by model configuration.", "Shading indicates the middle 50\\% (dark), middle 95\\% (medium),", "and entire range (light) of ensemble simulations", "(NPP/LAI results are clipped to (-5, 40) for readability,", "as a few isolated ensemble members produced extreme values at specific", "time points).", "Effective number of PFTs was calculated as the inverse", "of the Simpson diversity index with proportions expressed as", "fraction of aboveground biomass.", "Each colored line represents the same set of parameters across all configurations,", "and corresponds to the letters in Figure 5.", "Black dot with error bars is the observed mean and min/max value." # TODO REF ) fig("summary-ts-plot.png")
r pgbreak()
### Figure 9 lai_pft_cap <- paste( htmlfig(), "ED-2.2 PFT-level predictions of leaf area index (LAI) by model configuration and parameterization.", "Each row of plots is a simulation with the same set of parameters across model configurations.", "Letters correspond to the labels in Figure 7.", "Note that these are predictions of total leaf area index summed across all cohorts of the same PFT,", "meaning that higher LAI in one PFT does not necessarily indicate dominance, size, etc." ) fig("lai-pft.png")
r pgbreak()
### Figure 10 partial_variance_npp_cap <- paste( htmlfig(), "Contributions of individual parameters to predicted uncertainty in", "ED-2.2 predictions of average net primary productivity (NPP)", "over the period 1920-1950", "by model type.", "Each panel shows the top 8 parameter values by their fractional contribution", 'to parameter uncertainty ("partial variance").', "Input parameter distributions are shown in Figure 3.", "Note that x-axis scales vary across panels." ) fig("partial-variance-npp.png")
r pgbreak()
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.