Coupling earth system and integrated assessment models : the problem of steady state

Introduction Conclusions References


Introduction
Human activities are significantly altering biogeochemical cycles at the global scale, e.g. by appropriation of net primary production (Imhoff et al., 2004;Ito, 2011), modification of natural fire dynamics (Pechony and Shindell, 2010), and fossil fuel emissions raising atmospheric CO 2 levels (Le Queré et al., 2009) (LUC) exerts strong effects on the global carbon cycle (Bonan, 2008;Caspersen et al., 2000;Arora and Boer, 2010;Laganière et al., 2009), as well as direct biophysical effects on albedo and water vapor fluxes, that in turn have significant regional to global consequences (Brovkin et al., 2013;Jones et al., 2013b).As a result, different policy choices vis-à-vis LUC and carbon are projected to result in greatly different configurations of the future carbon cycle and climate system (Wise et al., 2009;Jones et al., 2013a).This raises a significant problem for global earth system models (ESMs), in which fully coupled climate models are used to draw inferences about Earth's past and future climate states and understand how changes to the radiative properties of Earth's atmosphere interact with its climate and biogeochemistry.Such models may incorporate static LUC inputs, but do not actively, or interactively, simulate policy options or economic forces, except for example in simple (in a policy sense) "thought experiments" (e.g., Bonan et al., 1992).This is a significant limitation given how strongly humans can perturb the earth system (Hurtt et al., 2002;Randerson et al., 2009).Conversely, integrated assessment models (IAMs) are used to examine the human components of the Earth system, including greenhouse gas emission sources, and drivers of landuse change.Their representation of the physical climate and earth system is simplistic, however, with little spatial resolution or process fidelity compared to an ESM (Meinshausen et al., 2011).These two modeling paradigms -ESMs with no economic or energy system modeling, and IAMs with only basic representations of natural processes -developed largely independently of each other, and their interactions have historically been limited.
ESMs and IAMs increasingly need each other's capabilities, however (van Vuuren et al., 2012;Houghton, 2013).One solution is to couple an ESM to an IAM, letting each model specialize in its specific domain while passing information on the natural and human systems, respectively, between them.This would provide a two-way coupling within a single integrated system, whereby economic decisions in the IAM translate directly into trace gas fluxes and land use changes in the ESM, and changes in the Introduction

Conclusions References
Tables Figures

Back Close
Full ESM climate would feed back onto crop yields, heating and cooling demands, energy production, etc. in the IAM.Successfully linking such complex, large models would permit integrated and unprecedented analyses of the interactions between economic change, climate policy, and the physical earth system, with fully coupled feedbacks between the economic and physical-science components (van Vuuren et al., 2012).
This paper describes the development and testing of such a mechanism linking an ESM (CESM, the Community Earth System Model) with an IAM (the Global Change Assessment Model, GCAM) (Fig. 1).The goals of the current study were to develop and test a robust but tractable coupling allowing GCAM to respond to changes in the CESM climate and biogeochemical cycles.Although we focus here on the CESM-to-GCAM coupling, this work is part of a larger effort to create a more general integrated Earth System Model (iESM) (Jones et al., 2013a) as described above.

Model descriptions
Both CESM's Community Land Model (CLM) and GCAM have been extensively described, and here we note only their aspects most relevant to this study (Gent et al., 2011).CLM version 4, used in this study, resulted from merging the biophysical framework of CLM v3.5 (Oleson et al., 2008) with the carbon and nitrogen dynamics of the biogeochemistry model Biome-BGC (Thornton et al., 2002;Running and Hunt, 1993).CLM incorporates biogeophysics, hydrologic cycle, biogeochemistry and dynamic vegetation (Bonan et al., 2002), and the dynamics of these model components have been extensively tested (Shi et al., 2011;Oleson et al., 2008;Lawrence et al., 2008;Mao et al., 2012a, b).Model vegetation is based on plant functional types (PFTs) occupying dynamic fractions of each grid cell (typically 0.25-2 • resolution), with each PFT (1 bare ground, 8 tree, 3 shrub, 3 grass, 1 crop) characterized by distinct physiological parameters (Oleson et al., 2010).The model's C and N cycles are closely coupled and include Introduction

Conclusions References
Tables Figures

Back Close
Full assimilation, plant growth and mortality, allocation, and subsurface cycling (Thornton et al., 2007); at any point in time, CLM tracks a wide suite of above-and belowground C pools resulting from the integrated effects of these and other (Kloster et al., 2010) processes.
The GCAM model, by contrast, is a dynamic recursive economic model driven by assumptions about population size and labor productivity that determine potential gross domestic product in each of 14 regions; these regions are further divided by GCAM's agriculture and land-use submodel into 18 agro-ecological zones or AEZs (Monfreda et al., 2009).GCAM originated as the energy-economic MiniCAM model (Edmonds and Reilly, 1983), and currently integrates energy, agriculture, forestry, and land markets with a simple terrestrial carbon cycle (Thomson et al., 2010;Wise et al., 2009).The model operates on a 5-year timestep, computing simultaneous market-clearing prices for all energy, agriculture, and land markets (Kim et al., 2006) under optional policy scenarios that can include carbon prices, emissions constraints, or capped limits on total radiative forcing (Calvin et al., 2009).Land cover, land use, agricultural and forestry production, and terrestrial carbon emissions are simulated endogenously, with initial above-and below-ground terrestrial carbon stocks (the only pools currently tracked) based on IPCC (2001) and Houghton (1999) data.Carbon emission and sequestration result from changes in land use between model simulation periods, with plant growth, decay, and soil C changes occurring at various rates depending on AEZ.
In the iESM architecture a third model, the Global Land Model or GLM (http: //gel.umd.edu/glm.php),currently downscales GCAM's land use decisions (made on agro-ecological zones at the regional level) onto CLM's global grid (Fig. 1).This step uses algorithms and assumptions described by Di Vittorio et al. (2014) and Lawrence et al. (2012), and is not detailed further here.

Linking the CLM and GCAM carbon cycles
The fundamental conceptual, as opposed to technical, problem in linking the CESM/CLM and GCAM models is that the former tracks real-time (actual) C pools and Introduction

Conclusions References
Tables Figures

Back Close
Full fluxes, while the latter bases its economic optimization on long-term (potential) C pools for large regions, and only computes changes -i.e., carbon emissions or uptake -due to LUC.Replacing GCAM's entire internal carbon cycle (and its reliance on potential carbon) may be possible in the long term, but in this study a looser coupling between the models was deemed more tractable, while also scientifically sufficient for the experiments described here.Such an approach transmits relative changes between the models while allowing baseline data, against which the models have been calibrated and tested, to differ.This means that when a CLM grid cell's carbon cycle changes, we need to (i) have a suitable proxy by which to change GCAM's steady-state carbon assumptions, as GCAM cannot currently use real-time carbon values, and (ii) distinguish LUC effects from climate and other (CO 2 , N deposition, etc.) effects, because only the latter should affect GCAM's decision-making.For example, if the carbon stock of a CLM forest changes from one time step to the next because of harvest, this should not affect GCAM's economic optimization -the forest will regrow to the same equilibrium state.If the same forest's carbon pool rises because of CO 2 fertilization, however, this information (i.e., there is more C sequestration potential available for this land use type) needs to be propagated to GCAM's assumptions about long-term pool potentials.Distinguishing these sources is thus critical (Gasser and Ciais, 2013), and these two issues are treated in turn below.

Proxies: a single-forcing test of CLM
Perhaps the most obvious proxy variable available to pass from CESM to GCAM is CLM's real-time carbon stock data, under the assumption that short-term stock changes will translate to longer-term storage changes.These data may be more vulnerable to LUC effects than carbon flux data, however, as fluxes typically recover much faster from disturbance than do the slower pools (Amiro et al., 2010;Goetz et al., 2012).
Environmental changes in C fluxes should in theory be related to steady-state C pools, even in the presence of landscape-scale disturbances (Hurtt et al., 2010), although the practical truth of this needed to be shown for CLM.Introduction

Conclusions References
Tables Figures

Back Close
Full We tested the feasibility of these proxies in two ways.First, we ran a series of single forcing factor experiments in CLM.The three factors used were atmospheric CO 2 , as alleviating the CO 2 constraints on leaf-level photosynthesis may cascade up to ecosystem carbon storage (Gedalof and Berg, 2010;Lenton and Huntingford, 2003); nitrogen deposition, a potentially strong constraint on the current and future global carbon cycle (Galloway et al., 2005;Norby et al., 2010); and LUC, which affects both immediate and long-term land-atmosphere interactions (Caspersen et al., 2000;Pongratz et al., 2009).
We performed four simulations to assess the relative contributions of these factors on output that could serve as potential proxy variables for GCAM: gross primary production, net primary production (NPP), heterotrophic respiration (HR), soil organic matter, vegetation carbon, and total ecosystem carbon.The CRUNCEP data used to drive these simulations is a combination of the CRU TS.2.1 0.5 • monthly 1901-2002 climatology (Mitchell and Jones, 2005) and the 2.5 • NCEP2 reanalysis data beginning in 1948 and available in near real time (Kanamitsu et al., 2002;Mao et al., 2012b).In allow the carbon stocks to approach the equilibrium state they would achieve in the absence of LUC.However, it is important to note that we did not disable the fire algorithms in CLM.Fire significantly influences model stocks and fluxes (Li et al., 2013), and thus rather than converging to a single steady-state carbon stock, PFTs influenced by fire converged to a quasi-equilibrium characterized by periodic carbon losses due to fire followed by periods of recovery.We examined the degree to which (i) NPP in the first 5 years of simulation S5 predicted total vegetation carbon in the final 5 years, and (ii) the change in NPP resulting from an altered climate state (S6 minus S5) predicted the relative change in C pools over the final years of the two simulations.
The overall control run-our base reference simulation-was based on the RCP4.5 stabilization pathway (Thomson et al., 2011), one of a series of Representative Concentration Pathways (Moss et al., 2010).The RCP4.5 scenario stabilizes global radiative forcing at 4.5 W m −2 in 2100, factoring in the effects of greenhouse gases, short-lived species, and LUC emissions.This is a cost-minimizing pathway, and reaching this target drives significant changes in the energy system and land use, with forests expanding from their present-day extent.This scenario was described by Thomson et al. (2011).More generally, the scenarios used here are replications of the RCP4.5 using an updated release of GCAM version 3.0 (Wise and Calvin, 2011;Kyle et al., 2011) (http://www.globalchange.umd.edu/models/gcam/download/).The updated model incorporates an agriculture and land-use component with higher spatial resolution, features a shorter time-step, and has modest changes to the underlying model assumptions.

Distinguishing climate from land-use signals
It is important to distinguish LUC and environmental changes when modeling or tracking the terrestrial carbon cycle (Houghton, 2013).For the iESM, even a perfect proxy variable will be subject to both climate and LUC during a CESM run, both before the run starts (i.e., during spinup or initialization phases) as well as during an experiment.For example, a cell in which a new PFT is established immediately prior to an iESM run Introduction

Conclusions References
Tables Figures

Back Close
Full would have very low C stocks and NPP in the first timesteps; as its vegetation regrows, the cell would appear, to GCAM, to be undergoing enormous productivity increases.Conversely, significant expansion of a PFT (e.g., agriculture reverting to forest) during the iESM run might appear to have drastically lowered productivity, leading GCAM to redirect resources away from that region.In both cases, we need to exclude these cells from the final numeric scalars (i.e., the proxy variables signaling how much GCAM should adjust its assumptions of steady-state C) computation for two reasons: first, they will bias the computation of the scalars; second, GCAM will, if it sees out-of-line values, potentially pour more resources into those cells, leading to a runaway feedback.
(A negative feedback is also possible; both cases occur because the changed productivity alters the relative profitability of the different land uses, and profit maximization is the fundamental decision-making criterion in GCAM.)To distinguish LUC from climate signals, we assumed that climate changes will have a broad spatial distribution, either global or regional, while LUC will affect relatively small groups of cells in any particular timestep; this obviously may not hold in particular regions and points in time (Arora and Boer, 2010), but should be broadly true across the millions of data points (∼ 10 5 grid cells × PFT combinations) being output by CESM.Thus a statistical outlier test, comparing how much any particular cell's carbon cycle has changed relative to the start of the run, should be able to exclude cells whose GCAM scalars, i.e. the inferred change in long-term carbon density, fell significantly outside of the norm.To do so we used a method based on median absolute deviation (Davies and Gather, 1993), a robust (insensitive to outliers) measure of central tendency.The scalars were then mapped from CLM'S PFTs and grid cells to GCAM's land-cover types and AEZ regions, weighted by PFT area, land area in each grid cell, and cell area in the AEZ.This technique depends on the overall population mean not being overly perturbed, and thus will not work in extreme scenarios of mass deforestation (e.g., Bonan et al., 1992).An important question is how soon, under increasing amounts of LUC, bias (i.e., LUC effects masquerading as climate change to GCAM) will be introduced into Introduction

Conclusions References
Tables Figures

Back Close
Full the iESM model system.We used a Monte Carlo simulation, written in the statistical package R 2.15.1 (R Development Core Team, 2012), to examine how robust this outlier exclusion method would be to different levels of LUC, and what if any bias it might introduce to the GCAM carbon density values.For this exercise, 10 000 cells were simulated in which a constant +10 % climate-change effect on the scalars was presumed to be occurring (i.e., a combination of CO 2 , climate, etc., that was causing a 10 % rise in productivity and/or C density) (Jain and Yang, 2005).A LUC effect, ranging from −500 % to +500 % and affecting from 5 % to 95 % of the cells, was then additionally applied.The outlier exclusion test defined above was then calculated on the cells, and the resulting scalars (i.e.inferred climate signal) compared to the original known climate signal to quantify the difference how much error would be introduced into iESM under such circumstances.

Single-forcing tests: identifying the best proxy variables
Clear differences emerged between the potential proxy variables tested in CLM in response to three different forcing factors (Fig. 2).Most notably, carbon stocks were much more sensitive to LUC than were carbon fluxes.This result matches both theory (Odum, 1969) and a wide variety of field studies (Amiro et al., 2010;Goetz et al., 2012): carbon stocks are by their nature integrative and accumulate relatively slowly compared to their respective fluxes.Carbon fluxes are thus more robust to LUC, and because we need to distinguish LUC from climate and other environmental drivers that influence productivity (including CO 2 and nitrogen deposition, shown in Fig. 2), we used the two primary fluxes determining carbon balance (net primary production and heterotrophic respiration, NPP and HR) as proxy variables linking CESM to GCAM for the next set of experiments.Specifically, CLM NPP changes were used to scale, on a relative basis, Introduction

Conclusions References
Tables Figures

Back Close
Full GCAM's aboveground potential carbon stocks, while a combination of NPP and HR provided a relative scaling for GCAM's belowground carbon.
A second, related problem arising from the use of carbon stocks as proxy variables can be seen in Fig. 3.In this case a test coupling between CLM and GCAM, using carbon stocks to pass climate-change information, produced sharp and unrealistic changes from the RCP4.5 control run.(This occurred even when running the outlierexclusion protocol described above.)Global LUC emissions climbed throughout the 21st century in a departure from the RCP4.5 control, because a few CLM grid cells, located in GCAM's "Middle East" region, were subject to LUC at the end of CLM's spinup phase.As a result, their C stocks (and GCAM's estimation of their long-term potential C) increased rapidly in the early years of the model run, leading GCAM to pour more resources into these cells (because these cells' productivity appeared extraordinarily high, as described in the Methods).Increasing the area of newly planted bioenergy crops created an even stronger signal of rapidly increasing carbon stocks, exacerbating the original problem and causing GCAM to put even more resources into the region.
By the end of the century, GCAM was mistakenly growing a huge percentage of the world's bioenergy crops in the region, on a very small area of land (Fig. 3).

Correlation between NPP and equilibrium pools in CLM
Simulations S5 and S6 provided insight into the relationship between NPP and equilibrium C pools within CLM.NPP at the beginning of the S5 simulation was a good predictor of the equilibrium pools values at the end of the simulation (Fig. 4), although the slope of this relationship varied for different PFTs.It was also apparent that this relationship breaks down at very low NPP values for some PFTs.This result is consistent with ecological theory and observations, as freshly disturbed ecosystems require a period of initial growth before NPP stabilizes.These very low NPP values were reliably excluded by the outlier exclusion method discussed above and tested below.
We also found that the change in NPP resulting from an altered pattern of climate (comparing simulations S5 and S6) was a relatively good predictor of the subsequent Introduction

Conclusions References
Tables Figures

Back Close
Full change in equilibrium C stocks.Table 1 shows the slopes of the linear relationships between the change in initial NPP (simulation S6 minus S5) and change in equilibrium C for each PFT in CLM.The initial change in NPP was able to explain 20-92 % of the variance in the C pool change over the 21st century simulation.In general, NPP was a better predictor for relatively high-carbon forest ecosystems, as compared to grasses, shrubs, and crops.
To determine whether fire dynamics were responsible for some of the unexplained variance in equilibrium C pools, we performed the same analysis a second time, excluding PFT-gridcell combinations in which the cumulative carbon loss from fire over the 150 year S5 simulation exceeded 800 g C m −2 .This led to moderate improvements in the R 2 values in all PFTs except the two broadleaf evergreen PFTs, and moderate increases in the regression slopes, indicating that fire-influenced regions tend to have lower C values than others.This is consistent with both observations and CLM's general fire characteristics (Li et al., 2013).This finding -that the relationship between NPP and equilibrium pools is improved by excluding fire-influenced regions -suggests that fire dynamics and fire regime changes in response to climate change are important to account for when constructing simple proxies that can predict changes in future terrestrial carbon stocks based on evolving climatic and ecological conditions.

Distinguishing LUC from climate
The initial experiments thus established the best available variables to loosely couple CESM and GCAM.But how well could the coupling -specifically, statistically excluding CLM grid cells whose carbon fluxes were changing "too fast" -separate LUC and climate signals?The Monte Carlo experiment results (Fig. 5) suggested that as long as fewer than ∼ 25 % of the simulation cells were perturbed, the error (between the known climate signal and that inferred by the outlier test) remained small (< 10 %).Even when larger numbers of cells were perturbed, the LUC effect had to be quite large to exceed this level.Because the iESM's outlier test is applied to the global population, and not sub-regions, this implies that only under extreme scenarios will this mechanism start Introduction

Conclusions References
Tables Figures

Back Close
Full to introduce substantial error.(In test iESM runs attempting to reproduce RCP 4.5, 4-8 % of the global grid cells were excluded -i.e., failed the outlier test -at each timestep.)Such scenarios can be a useful tool, however (Bonan, 2008;Nobre et al., 1991;Thomson et al., 2010), and a full integration of the iESM models' carbon cycles is clearly desirable in the future.

Conclusions
Here we have implemented and tested a coupling mechanism between the carbon cycles of an earth system model (CLM) and an integrated assessment (GCAM) model.CLM's net primary production and heterotrophic respiration outputs were found to be the most robust proxy variables by which to manipulate GCAM's assumptions of longterm ecosystem steady state carbon, with short-term forest NPP shifts strongly correlated with long-term biomass changes in particular.By assuming the carbon cycle effects of land-use change are short-term and spatially limited relative to widely distributed climate effects, we were able to distinguish these effects successfully in the model coupling, passing only the latter to GCAM.Increasingly extreme LUC scenarios will eventually break down this mechanism, however.By allowing climate effects from a full earth system model to modulate, in real time, the economic and policy decisions of an integrated assessment model, this work provides a foundation for further development of the iESM project linking these models in a robust and flexible framework.More generally, by allowing climate effects from a full earth system model to dynamically modulate the economic and policy decisions of an integrated assessment model, this work provides a foundation for linking these models in a robust and flexible framework capable of examining two-way interactions between human and earth system processes.Introduction

Conclusions References
Tables Figures

Back Close
Full  Full . In addition, land-use change Figures Back Close Full Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | simulation S1 (the control), we used 1901-1920 climate drivers for the entire period 1850-2010, and kept atmospheric CO 2 concentration, nitrogen deposition, and land use boundary conditions constant at their 1850 values.In simulations S2-S4, we used the same looped climate, and varied one of the three factors in each.The effect of each individual factor was then calculated by subtracting S1 from simulations S2, S3 and S4.Second, we examined how well NPP in particular was related to equilibrium C stocks in CLM.This involved two offline experiments in which land cover was fixed at 2000 levels, with a repeating 5-year climate drawn either from the beginning (2005-2009, simulation S5) or end (2090-2094, simulation S6) of an RCP4.5 coupled simulation performed at NCAR as part of the 5th Coupled Model Intercomparison Project (Taylor et al., 2012).The state of the terrestrial carbon system at the beginning of these simulations reflected the disturbance and climate histories of the 20th century, with a spectrum of different non-equilibrium states for various PFTs in different grid cells.We ran each simulation for 150 model years with no additional land-use change in order to Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | entific Computing Center, which is supported by the Office of Science of the US Department of Energy under Contract DE-AC02-05CH11231.The CESM project is sup-ported by the National Science Foundation and the Office of Science (Biological and Environmental Research) of the US Department of Energy.We thank S. Smith for his thoughtful comments on an early draftDiscussion Paper | Discussion Paper | Discussion Paper |

Fig. 5 .
Fig. 5. Monte Carlo simulation examining how well an outlier test can distinguish between climate and land use change (LUC) signals when passing data from CLM to GCAM.Contour lines (every 25 %) show error between the inferred climate change signal and known (artificial) signal as increasing numbers of cells (y axis) are perturbed by LUC with increasing intensity (x axis).

Table 1 .
Slope (yr), adjusted R 2 value, and number of grid cells for the relationship between change in NPP in response to a climate change signal and resulting change in equilibrium biomass.Excluding PFTs whose cumulative carbon loss from fires exceeds 8 Mg C ha −1 over 150 years generally improved the R 2 values and increased the slopes (data not shown).