Contribution of permafrost soils to the global carbon budget

For this study, the dynamic global vegetation model LPJmL (Sitch et al 2003, Bondeau et al 2007) has been extended by a new permafrost module. This permafrost implementation is based on a detailed energy balance, explicit freezing/thawing mechanics and an improved soil hydrology model, which goes beyond the previous implementations of Beer et al (2007) and Wania et al (2009). LPJmL explicitly simulates key ecosystem processes such as photosynthesis (Farquhar et al 1980, Collatz et al 1991, 1992, Haxeltine and Prentice 1996a, 1996b), plant and soil respiration, carbon allocation, evapotranspiration and phenology of 9 plant functional types (PFTs) representing natural vegetation at the level of biomes (Sitch et al 2003), and of 11 crop functional types (Bondeau et al 2007). The water and carbon cycles are fully coupled (Gerten et al 2004).

To better simulate the carbon dynamic in the permafrost zone, which is affected by vegetation growth and soil carbon decomposition, a number of model modifications have been implemented. The soil carbon decomposition rate is mainly determined by soil hydrology and soil temperature. The introduction of a newly implemented permafrost module and a new hydrology scheme has led to an improved representation of the decomposition process. Figure S2 (available at stacks.iop.org/ERL/8/014026/mmedia) shows that the model is able to simulate net ecosystem exchange (NEE) at a number of boreal fluxnet sites, which is calculated from the difference of heterotrophic respiration and net primary production (NPP). To account for different decomposition rates in the different soil layers, which is of particular importance in permafrost zones, a vertical distribution of soil carbon within the soil column following Jobbagy and Jackson (2000) has been introduced to the model. The total soil carbon amount SoilC_total is estimated from the mean annual decomposition rate k_mean(i) and the mean litter input into the soil as in Sitch et al (2003), but is distributed to all root layers separately (equation (1)). The envisaged vertical soil distribution SoilC_(i):

$$\begin{eqnarray} &&{\mathrm{SoilC}}_{(i)}=\sum _{p=1}^{n p f t}{d}_{(i)}^{{K}_{(p)}}{\mathrm{SoilC}}_{\mathrm{total}} \end{eqnarray} \tag{ 1 }$$

is estimated after a carbon equilibrium phase of 2310 years. d is the relative share of the layer i in the entire soil bucket. The slope parameters K characterize the relative rate of decrease of soil carbon content with depth for each specific PFT (p). Values are shown in table S2 (available at stacks.iop.org/ERL/8/014026/mmedia). The vertical soil distribution shown in table 4 of Jobbagy and Jackson (2000) is characterized by the soil carbon share of the respective layer and PFT socfr_(i,p):

$$\begin{eqnarray} &&{\mathrm{socfr}}_{(i,p)}=1{0}^{({K}_{(p)}\log 1 0({d}_{(i)}))}. \end{eqnarray} \tag{ 2 }$$

The mean decomposition rate for each PFT k_meanPFT(p) can be derived from the mean annual decomposition rate k_mean(i) of the spinup years as a layer weighted value derived from equation (2):

$$\begin{eqnarray} {k}_{{{\mathrm{mean}}_{\mathrm{PFT}}}_{(p)}}&=&\sum _{i=1}^{n\mathrm{layer}}({k}_{{\mathrm{ mean}}_{(i)}}{\mathrm{socfr}}_{(i,p)}) \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} {C}_{{\mathrm{shift}}_{(i,p)}}&=&\frac{{\mathrm{socfr}}_{(i,p)}{k}_{{\mathrm{mean}}_{(i)}}}{{k}_{{{\mathrm{mean}}_{\mathrm{PFT}}}_{(p)}}}. \end{eqnarray} \tag{ 4 }$$

The annual carbon shift rates C_shift(i,p) describe the organic matter input from the different PFTs into the respective layer due to cryoturbation and bioturbation and are designed for global applications.

The thermal soil characteristics are described by a heat conductance energy balance model, accounting for latent heat from freezing and thawing of soil water:

$$\begin{eqnarray} &&\frac{\partial {T}_{(x,t)}}{\partial t}=\alpha \frac{{\partial }^{2}T}{\partial {x}^{2}} \end{eqnarray} \tag{ 5 }$$

where α is thermal diffusivity, and T soil temperature at position x and time t.

Soil humidity is calculated dynamically in response to thawing and freezing, percolation of water above field capacity and transpiration of plant available water in the upper 3 m of the soil column. The first 20 cm are influenced additionally by soil evaporation. Freezing depth is calculated assuming that the fraction of frozen water is congruent with the frozen soil bucket. All newly implemented processes, which are described and evaluated in detail in the supporting information (SI available at stacks.iop.org/ERL/8/014026/mmedia), enable the dynamic consideration of soil carbon dynamics in regions affected by permafrost. Soil data describing the thermal and hydraulic characteristics are taken from the Harmonized World Soil Database (version 1.2) (2012). The original data set was first aggregated to 0.5° resolution and then classified according to the USDA soil texture classification (http://edis.ifas.ufl.edu/ss169). Hydraulic soil parameters are derived from Cosby et al (1984) and thermal parameters from Lawrence and Slater (2008) (see table S1 available at stacks.iop.org/ERL/8/014026/mmedia).

The model is capable of reproducing local carbon fluxes and water flux dynamics, as demonstrated in a comparison with eddy-flux tower measurements (see supplementary figures S2 and S3 available at stacks.iop.org/ERL/8/014026/mmedia). For Ameriflux and Euroflux ORNL DAAC (2011) sites of the northern latitudes above 50°N the model is able to reproduce the size and most of the variability in observed carbon and water fluxes. The observed runoff spring peak (resulting from snow melt and precipitation on still frozen ground) and interannual runoff variability in the permafrost affected river catchments is also reproduced well by the model (figure S4 available at stacks.iop.org/ERL/8/014026/mmedia). In some cases the timing is not met, which is a result of the constant value of flow velocity of 1 m s⁻¹ used (Rost et al 2008). The quality similarity index (QSI) shows the similarity of the time series (figure S4 available at stacks.iop.org/ERL/8/014026/mmedia), especially if a time shift occurs between observation and simulation (Jachner et al 2007). The resulting soil carbon distribution calculated by LPJmL is shown in figure S5 (available at stacks.iop.org/ERL/8/014026/mmedia).

The Climate Research Unit's (CRU) time series (TS) 3.1 data for temperature and cloud cover (Mitchell and Jones 2005) were used for model spinup and the subsequent 20th century transient run. For precipitation the Global Precipitation Climatology Centre's (GPCC) gridded precipitation data (version 5) (Rudolf et al 2010) was used. As the land mask for the two data sets differ, the GPCC data set was adjusted to the CRU land mask by interpolating missing cells (Heinke et al 2012). The number of wet days per month was created synthetically based on precipitation data using the approach by New et al (2000).

For the 21st century simulations climate scenarios were generated using a pattern-scaling approach (Heinke et al 2012). Climate patterns describing the GCM-specific changes of local temperature, precipitation and cloud cover in response to changes in global mean temperature (GMT) were extracted from the 19 GCMs which supplied data for the SRES scenarios to the World Climate Research Program's Climate Model Intercomparison Project (WCRP CMIP3). These climate patterns were scaled using GMT trajectories to create transient climate anomalies over the 21st century. GMT and corresponding greenhouse gas trajectories for temperature increases of 1.5–5.0 K above pre-industrial level by 2100 in steps of 0.5 K were generated using the reduced complexity climate model MAGICC6 (Meinshausen et al 2011). Transient climate anomalies were applied to a baseline climatology based on observations that deliver interannual variability. Wet days were derived synthetically from precipitation data for the historical period following the approach of New et al (2000). For a detailed description see (Heinke et al 2012).

The model runs on 0.5° × 0.5° spatial resolution using monthly input of the above described data and with global annual atmospheric CO₂ concentration from Mauna Loa station (NOAA/ESRL, www.esrl.noaa.gov/gmd/ccgg/trends/) for the historical period and with the derived CO₂ concentration from MAGICC6 for future scenarios. The model runs on a daily time step by linear interpolation of monthly data and for daily precipitation with a generator based on wet days and precipitation (Gerten et al 2004). The carbon equilibrium was reached after 5000 years of simulation recycling the first 30 years of the 20th century with the observed climatology as described above. Then additional 390 years were used to simulate the effects of historic land use change (Fader et al 2010) on the carbon budget until 1900, recycling again climate data of the years 1901–30, followed by the transient run with observed climate and land use data until 2009. As the land use data start in 1700, for the spinup time 1511–700 the land use distribution of 1700 was taken. Areas under cultivation for crop land and grazing land were included (though they are negligible in the permafrost zone), as they are required for a correct assessment of the present global carbon budget and the contribution of the permafrost zone to it. Land use was assumed to remain constant for the scenario runs.

All 152 climate scenarios (19 GCM patterns, 8 GMT trajectories) were started from the same conditions in 2010, reaching the respective GMT increase in 2100 (average of 2086–115). These runs (referred to as 'PF') were extended until 2615 with constant temperature levels and a constant atmospheric CO₂ content from the year 2100 in order to analyse equilibrium times and soil carbon stocks of long-term stabilization. A second experiment, ('PF-MH'), was conducted with the same procedure as the PF run, but assuming microbial heating of 40 MJ kgC⁻¹ as suggested by Khvorostyanov et al (2008b) and used by Koven et al (2011). A third experiment ('PF-Y') additionally accounts for deeper soil carbon deposits in soil layers below 3 m (the so-called yedoma deposits in Eastern Siberia). These were assumed to contain the same amount of carbon as the slow carbon pool of each grid cell in 2010 to accommodate the spatial distribution for the scenario period. The decomposition rate of the deep soil organic matter was parametrized with a rate of 0.025 yr⁻¹ following Dutta et al (2006) and Zimov et al (2006a).

An additional experiment was conducted where change in biomass input from changed plant growth and vegetation distribution was neglected. This was implemented by running the model with static biomass input into the litter pool, allocating the average litter input of the historical period (1980–2009) to each cell and year.

The model results are in good agreement with independent studies of the current permafrost distribution and its active-layer thickness (figure 1). The study area, i.e. the permafrost zone, is defined to be north of 50°N and bounded by a maximum thaw depth of less than 13 m. For soil carbon stored in the upper 3 m of the permafrost zone the reported range spans from 400 Pg in a model study by Koven et al (2011) to 1024 Pg given by an up-scaling approach of different soil databases (Tarnocai et al 2009). This study estimates 952 Pg carbon in 2009, which is in the upper range, but more recent studies tend to give higher estimates (Schuur et al 2008). Furthermore the deep soil carbon yedoma deposits include additional 407 PgC in the deeper layer (Tarnocai et al 2009). The LPJmL model estimates additional 343 PgC for the PF-Y experiment, in good agreement with estimations by Zimov et al (2006a, 2006b), and Schuur et al (2008) of carbon contained in deep loess sediment accumulations below 3 m. Taking additional microbial heating into account (PF-MH), the amount of soil carbon in the pre-historic equilibrium state is quite similar to that of the PF run (carbon equilibrium was calculated as described above). Contrary to expectations, the soil carbon pool is slightly higher in PF-MH than without microbial heating, about 10 PgC in the overall permafrost zone. This is due to the thicker active layer leading to more available water for plant transpiration, and correspondingly higher plant productivity. We here focus on the experiment without additional heating as Koven et al (2011) have pointed out that this process shows a strong threshold behaviour and a large uncertainty due to a lack of quantitative observations.

Zoom In Zoom Out Reset image size

The global land biosphere uptake modelled by LPJmL largely explains the interannual variability of the global carbon sink (figure 2). By subtracting emissions from fossil fuel combustion and cement production from the annual growth rate of atmospheric CO₂ concentration the strength of the apparent global sink can be quantified. The inferred residual terrestrial sink including land use change fluxes and ocean fluxes as suggested by Le Quéré et al (2009) and the terrestrial sink modelled by LPJmL (including land use change emissions) are of similar magnitude and are highly correlated. The years 1991 and 1992 are the only years which cannot be explained by here simulated biosphere dynamics. Only a small increase of the atmospheric CO₂ concentration is measured in these years. Possible explanations are that the effects of the Pinatubo eruption are not included in both models or that the emission data are inaccurate due to underestimation of the collapse of industrial activity in the former eastern bloc countries especially in the year 1991 (Peters et al 2012). However, in the years 2000–9 both estimates of the carbon sink are in a good agreement. We estimate a loss of 7.5 Pg soil carbon globally during this period. The overall sink is largely determined by the increase of vegetation carbon which is estimated to be about 16 PgC in the years 2000–9.

Zoom In Zoom Out Reset image size

More than half of the global sink is located in the permafrost zone (table 1) where our simulations show large carbon stocks. A loss of this sink would strongly affect atmospheric CO₂ concentration as less of the anthropogenic emissions would be taken up by the terrestrial biosphere. The magnitude of the land–atmosphere flux, the balance of soil respiration, fire carbon emissions and NPP, in the permafrost zone has been dominated by the increase in vegetation carbon (table 2) since the 1980s. Soil carbon has been increasing as well, but the uptake rate has been weakening. In the years 2000–2009 the soil carbon pool was almost static (see table 2) but has not yet turned into a carbon source.

Average temperatures in the study area have increased by 1.5 K since the beginning of the 20th century, which is more than the GMT increase. All climate scenarios used in this study indicate that this phenomenon of so-called polar amplification will continue during this century (table S2 available at stacks.iop.org/ERL/8/014026/mmedia). We find that for most GMT levels simulated carbon storage in permafrost soils decreases until 2100 (figure 3 and table 3). The GMT-dependent changes in soil carbon can be explained by the following mechanism. Due to increased plant productivity in response to CO₂ fertilization, longer growing periods and migration of trees into the permafrost zone, the organic matter input into the soil increases. Modelled vegetation carbon pools in the study area increase by 9–70 PgC by 2100. This is followed by an accumulation of soil carbon of 1.7–8.3 Pg in the first half of the century. After 2060, however, soil carbon pools decrease because the effect of rising temperatures on soil decomposition becomes stronger than the effect on litter input. Carbon uptake and release equalized only in the lowest temperature change scenario of 1.5 K. The slow carbon pool declines for all temperature levels, which indicates that the enhanced decomposition cannot be compensated for by increased litter input. The probability that permafrost soils turn into a carbon source is higher than 66% over all climate scenarios (including all GMT increases) and rises to 96% in both PF-Y and PF-MH experiments. The simulation of microbial heating leads to a positive feedback that increases thawing depth and accelerates the decomposition of soil organic matter. Furthermore the decomposition of the deep soil carbon deposits will lead to a much higher carbon release of more than 9 PgC until 2100 including 90% of all climate scenarios, and the median response of the 5.0 K GMT level is estimated to be more than 22 PgC. If vegetation dynamics is not taking into account the soil carbon release will be much faster (figure S6 available at stacks.iop.org/ERL/8/014026/mmedia). All scenarios start to emit carbon from the beginning of the scenario time and in 2100 at least 27 PgC is lost and carbon decrease for 5.0 K GMT increase ranges from 75 PgC to 120 PgC (figure S6 available at stacks.iop.org/ERL/8/014026/mmedia).

Zoom In Zoom Out Reset image size

Additional carbon input into the soils is simulated to occur predominately in the far northern regions (figure 4) where the tundra vegetation mostly represented as C3 grasses becomes more productive. North of 65°N the root zone does not thaw completely even under 5.0 K GMT increase, hence a large amount of carbon is still frozen and not released into the atmosphere until 2100. In the study area south of 60°N mean soil temperatures exceed 0 °C under all levels of GMT increase and large areas are no longer dominated by permafrost. The loss of permafrost area in 2100 at 2 K is estimated to be about 12% of the current area for 90% of all climate scenarios and twice as much for the 5 K scenarios is very likely. As a consequence, NPP increases due to increased rooting depth and higher water availability and due to an increased photosynthesis rate. Grasslands benefit most from this NPP enhancement, but this only causes a slight increase of vegetation carbon (figure 4). Trees are not able to benefit, rather a decline of tree cover is projected south of 60°N, caused by higher atmospheric water demand which cannot be fulfilled owing to the limited rooting depth of trees. Between 60°N and 75°N trees do benefit from the temperature increase, and the CO₂ fertilization effect and resulting NPP enhancement, leading to an increase of tree cover that at the beginning of the 21st century could only be found south of 60°N. Consequently, tree cover is rather homogeneous from 50°N to 70°N under the scenarios of high GMT increase.

Zoom In Zoom Out Reset image size

In our simulations, the permafrost zone loses its capacity to sequester carbon under almost all considered levels of global warming (figure 5). The current sink will be amplified until 2050 but already around 2080 it will diminish or even turn into a source in 2100. Only for 3.0 K GMT increase the permafrost zone still constitutes a small sink. This is due to the increased NPP and reduced soil carbon loss at this GMT level. The lower GMT increases induce a smaller increase of NPP, for 1.5 K even a decrease, and a lower soil carbon release compared to the higher GMT levels. Both mechanisms lead to a net carbon loss of terrestrial carbon in 2100. If the soil carbon input does not change, the land–atmosphere fluxes already show a source at the beginning of the 21th century (figure S7 available at stacks.iop.org/ERL/8/014026/mmedia). The carbon source strength in 2100 reaches at least 0.25 Pg yr⁻¹ and for the highest temperature scenario the carbon source rises up to 4.4 Pg yr⁻¹.

Zoom In Zoom Out Reset image size

The simulation experiments for a stabilized warmer temperature show that the benefit of enhanced plant growth will not compensate for the soil carbon loss even over a longer period and the soil carbon equilibrium will be achieved at lower carbon contents for all GMT levels (figure 6). We also find that soil respiration increases for all levels of GMT increase. An equilibrium occurs much faster under low levels of warming and in the 5.0 K increase scenario soil carbon pools still decrease at the end of the simulation period after additional 500 simulation years. Consequently, all levels of GMT increase will result in a net carbon loss (figure 6), due to the delay of the thermal response in permafrost soils. In response to further thawing in the long-term, the deep soil carbon deposits (more than 3 m) are subjected to soil carbon release. This means that the response of permafrost soils to temperature increase will cause a long-term signal.

Zoom In Zoom Out Reset image size

The different GMT levels also have distinct effects on vegetation and soil carbon pools (figure 7). Soil carbon storage decreases already at 1.5 K warming, although plant productivity and consequently carbon input into the soil increase at the same time. The permafrost soil carbon pool is diminished by 14 Pg at 1.5 K GMT increase and by 178 Pg at 5.0 K GMT increase for the stabilization experiment. The sensitivity of permafrost soil carbon to global warming increases non-linearly with higher temperature. As soon as vegetation becomes less productive at levels between 4.5 and 5.0 K temperature increase, vegetation carbon pools decline and soil carbon pools decline much faster than under lower temperature increases (figure 7).

Zoom In Zoom Out Reset image size