Model estimates of climate controls on pan-Arctic wetland methane emissions

Introduction Conclusions References


Introduction
Methane (CH 4 ) is an important greenhouse gas, with a greenhouse warming potential about 25 times that of CO 2 (IPCC, 2013).Globally, wetlands are the largest natural CH 4 source (Fung et al., 1991;Hein et al., 1997;IPCC, 2013).The strong sensitivity of wetland CH 4 emissions to ambient soil conditions has led to concerns about possible feedbacks to climate change (Gedney et al., 2004;Eliseev et al., 2008).The northern high latitudes contain about one half of the world's wetlands (Lehner and Döll, 2004) and are experiencing more rapid climate change than elsewhere globally (Serreze et al., 2000;Diffenbaugh and Giorgi, 2012).The potential liberation of vast quantities of carbon from thawing permafrost provides additional impetus to efforts to understand the sensitivity of northern wetland CH 4 emissions to climate change (Schaefer et al., 2011;Koven et al., 2011).CH 4 emission rates in northern wetlands (which are predominantly peatlands) depend on a number of environmental and climate controls, including soil temperature, water table depth, labile carbon substrate, soil pH, oxidation state, nutrient concentrations, and vegetation composition (Saarnio et al., 1997;Christensen et al., 2003;Zhuang et al., 2004;Riley et al., 2011;Spahni et al., 2011;Glagolev et al., 2011;Lupascu et al., 2012;Levy et al., 2012;Olefeldt et al., 2013;Sabrekov et al., 2014).Many of these factors can interact and compete.For example, Bohn et al. (2007) showed via a process-based model that air temperature and precipitation exert competing influences on (a) water table depth, through winter snow accumulation, spring snowmelt, and summer precipitation and evapotranspiration, and (b) metabolic rates, through soil temperature, leading to trade-offs in their influ-Published by Copernicus Publications on behalf of the European Geosciences Union.

X. Chen et al.: Model estimates of climate controls
ences on emissions.Extreme (limiting) values of one factor can raise the sensitivity of emissions to that factor (Olefeldt et al., 2013).As a result of these interactions, different factors exert dominant controls at different sites (Olefeldt et al., 2013) or timescales (Sabrekov et al., 2014), hindering efforts to constrain model behaviors in the face of sparse observations (Melton et al., 2013).Therefore, isolating those conditions under which different factors dominate or limit the response of wetland methane emissions to climate change would benefit future field campaigns and modeling studies.
Previous attempts to characterize the sensitivities of northern wetland CH 4 emissions to environmental factors have included both data-driven (Bloom et al., 2010;Olefeldt et al., 2013) and process-based modeling (Bohn et al., 2007;Ringeval et al., 2010) approaches.Data-driven studies have the potential advantages of relative accuracy and simplicity but can have limited predictive power.For example, Olefeldt et al. (2013) found clear relationships between observed emissions from over 300 high-latitude sites and soil temperature, water table depth, and vegetation composition.However, while these relationships are a crucial step forward in our understanding, they must be embedded within a processbased model to estimate the aggregate response of northern wetland emissions to a given change in climate or characterize how these relationships may change with changing climate.Bloom et al. (2010) fit a regression model to observed atmospheric CH 4 concentrations from the Scanning Imaging Absorption Spectrometer for Atmospheric Chemistry (SCIAMACHY; Bovensmann et al., 1999) to observed surface temperatures from the National Center for Environmental Prediction/National Center for Atmospheric Research (NCEP/NCAR) weather analyses (Kalnay et al., 1996) and gravity anomalies from the Gravity Recovery and Climate Experiment satellite (GRACE; Tapley et al., 2004) and found that air temperature exerted the dominant control over highlatitude emissions.Unfortunately, the short (4 years) record length and the use of GRACE data as a proxy for near-surface wetland soil moisture suggest that these findings are highly uncertain and limited to the time span of the satellite data sets used.
Process-based studies potentially have greater predictive power, but their relative complexity may involve highly uncertain parameterizations.For example, Ringeval et al. (2010) found that variations in inundated area contributed 30 % to the interannual variability in CH 4 emissions over the latitudes north of 50 • N.However, despite the strong emissions observed at non-inundated peatlands throughout the high latitudes (e.g., Saarnio et al., 1997;Panikov and Dedysh, 2000;Friborg et al., 2003;Glagolev et al., 2011), they only considered emissions from inundated wetlands, thus potentially inflating the contribution attributed to inundation.Bohn et al. (2007) accounted for non-inundated emissions, but their study was restricted to a small area in West Siberia.Numerous other process-based studies (using both forward and inverse models) have investigated the response of northern CH 4 emissions to historical or future climate variations (e.g., Chen and Prinn, 2006;Bousquet et al., 2011;Riley et al., 2011;Spahni et al., 2011;Bohn et al., 2013Bohn et al., , 2015;;Zhu et al., 2014), but none have attempted to characterize the sensitivities of emissions to climate factors as a function of geographic location, wetland type, or climate conditions.CH 4 emissions are not the only biogeochemical process for which environmental controls have been investigated.Nemani et al. (2003) found that annual net primary productivity (NPP) is limited by temperature and radiation at high latitudes but by moisture-related factors at lower latitudes.Teuling et al. (2009) and Seneviratne et al. (2010) investigated global climate controls on annual evapotranspiration (ET), and found that temperature is the dominant control over northern Eurasia, while precipitation is the dominant control at mid-to low latitudes and in northern Canada.These datadriven studies all produced maps of the regions in which various climate factors dominate the flux in question.Such maps are useful in understanding how climate factors interact, which processes are most important, and how these fluxes might evolve under future climate change, particularly in cases where observations are sparse (as is true for CH 4 emissions).
In this study, we use a process-based model to characterize the dominant climate drivers of northern high-latitude wetland CH 4 emissions and how they will change with changing climate.We address three questions: 1. What have been the aggregate long-term CH 4 emissions from the pan-Arctic wetland area over the last 50 years, and how have they changed?
2. What have been the dominant factor(s) controlling changes in the space-time variability of CH 4 emissions over that time period?
3. How will these conditions be affected by a changing climate over the remainder of the 21st century?
To investigate these questions, we use an enhanced version of the Variable Infiltration Capacity (VIC) large-scale hydrology model (Liang et al., 1994;Bohn et al., 2013) and the wetland CH 4 emissions model of Walter and Heimann (2000).In answering questions (2) and (3), we develop (a) maps of the sensitivities of simulated pan-Arctic wetland CH 4 emissions to various environmental factors, (b) maps of correlations between these factors and CH 4 emissions, and (c) empirical estimates of how these sensitivities and correlations depend on climate.These sensitivity maps and climate dependencies provide a basis for projecting future emissions in the region, which we then compare with our VIC model projections to evaluate their ability to capture the effects of underlying processes.

Spatial domain
Our study domain is the global land area north of 45 • N (Fig. 1a) with slight modifications.Because this region contains all the river basins that drain into the Arctic Ocean, we will refer to it as the "pan-Arctic" hereafter, as in Slater et al. (2007).Our domain boundaries are as in the TransCom project (Gurney et al., 2000), except that we exclude Greenland.We also include southern Russia and the permafrost part of Tibet.We divided the domain into 3775 100 km EASE (Equal-Area Scalable Earth) grid cells (Brodzik and Knowles, 2002).Our domain includes three major wetland areas (Lehner and Döll, 2004;Tarnocai et al., 2009;Fig. 1b): the West Siberian Lowland (WSL), which we define as the region from 55 to 75 • N and 60 to 90 • E; Scandinavia (55-75 • N and 15-45 • E); and the Hudson's Bay Lowland (HBL), which we define as the region from 45 to 60 • N and 75 to 100 • W.There are also many smaller wetlands distributed over the domain.The vast majority of the domain's wetlands are peatlands, which are reservoirs of organic carbon (Tarnocai et al., 2009), and have the potential to produce huge fluxes of carbon (CO 2 or CH 4 ) to the atmosphere.Forests cover about 23 % of the total land area of our study domain, as evidenced by the belt of high values of leaf area index (LAI) between about 55 and 65 • N (Myneni et al., 2002;Fig. 1c).
Our domain includes essentially the entire Northern Hemisphere permafrost land area, aside from a few high-altitude areas (Fig. 1d; see also Brown et al., 2014).Within the permafrost areas, deep soil temperatures are generally below 0 • C for successive years, which restricts biological methanogens.However, during summer, the active layer (seasonally thawed) provides a suitable environment for CH 4 production.

Model framework
We used a modified version of the Variable Infiltration Capacity (VIC) version 4.1.2(Liang et al., 1994;Bohn et al., 2013) that simulates carbon fluxes as well as the hydrologic processes represented in the standard version of the VIC model.The VIC model resolves the soil moisture and temperature profiles through a coupled water-energy balance scheme that accounts for cold-climate processes such as soil freeze-thaw and the insulating effects of organic soils.We provide here a brief description of the model features related to wetland process.The main enhancement in the version of VIC we used is a module for calculating the carbon inputs into the ecosystem, which is the substrate source of biogeochemical processes that produce CH 4 .Within each grid cell the model represents multiple land cover "tiles".This modified version of VIC also represents lakes and wetlands as described in Bohn et al. (2013).Each grid cell in the study domain is assumed to be composed of a lake-wetland tile and an upland portion (that may contain several different land cover tiles).The lake-wetland tile contains peatlands of fixed area, within which a time-varying portion may be seasonally inundated and which may contain a permanent lake.Peatlands, which are modeled as a mix of moss and shrubs, are allowed to emit CH 4 , subject to oxidation above the water table, but lake CH 4 emissions are set to 0. The water table depth within peatlands follows a distribution derived from assumed microtopography.Net primary productivity (NPP) within peatlands experiences inhibition when the water table is above the soil surface.More details of the lake-wetland continuum are included in Bohn et al. (2013).
Permanent lakes were prescribed using the Global Lakes and Wetlands data set (GLWD) of Lehner and Döll (2004).Wetland areas were taken in most cases from the union of wetland classes from the GLWD and wetland pixels from the MODIS plant functional type data set MCD12Q1 (Friedl et al., 2010).However, in regions where the GLWD delineated wetland classes as 25-50 and 50-100 % (occurring in Alaska and Canada), we defined wetlands as pixels with soil organic carbon content above 70 % from the Northern Circumpolar Soil Carbon Database (Tarnocai et al., 2009).Of the domain's 3775 cells, 2049 contain wetlands (lake-wetland fractions shown in Fig. 1b).
The enhanced VIC model is linked to the Walter and Heimann wetland CH 4 emissions model (Walter and Heimann, 2000), as described in Bohn et al. (2013) The combined VIC and CH 4 models were calibrated over West Siberia in Bohn et al. (2013), and we adopted the median parameter values from the distributions from that study (Table 1) for our primary simulations.In Bohn et al. (2013), two parameter sets were optimized for the West Siberian Lowland: "south" (primarily within the forest belt, or taiga) and "north" (primarily tundra).These parameter sets only corresponded to broad geographic regions, rather than to specific types of wetlands such as bogs or fens.To extend these parameter sets across our entire domain, we assigned the "south" parameter set to grid cells with July LAI higher than 4 and the "north" parameter set to all other grid cells.
LAI data were taken from the MODIS MCD15A2 data set (Myneni et al., 2002) for the period 2002-2010.We used the mean seasonal cycle for this period repeatedly for every year in our simulation period.Soil parameters were taken from Su et al. (2006).
The primary meteorological forcings used to drive the VIC include 3 h precipitation, air temperature, wind speed, downward shortwave and longwave radiation.These data were obtained from Sheffield et al. (2006) at 0.25 × 0.25 • spatial resolution, which we regridded to a 100 km EASE grid.Atmospheric CO 2 concentration data were taken from Bohn et al. (2013).

Simulations
Our historical simulation period was 1948-2006.Model spin-up consisted of two stages: (1) initialization of carbon pool storages and (2) a 50-year spin-up to stabilize moisture and carbon pools.We initialized soil carbon pools via an iterative procedure in which we identified the initial storage that would result in zero net change in carbon storage over the period 1948-1957.Then, to account for the pools' not yet having reached equilibrium with recent Holocene climate, we rescaled all three pool storages by the ratio of observed to simulated total carbon storage across West Siberia, using observations from Sheng et al. (2004).Then we ran the model for 50 years (5 × the decade 1948-1957) to stabilize its moisture and carbon storages.Starting from the model state at the end of this 50-year spin-up, we then performed simulations for 1948-2006.
To isolate the effects of various climate factors that drive the variability in CH 4 emissions, we performed five control experiments in which we removed trends (at each grid cell) in one or more variables (air temperature and longwave radiation; precipitation; air temperature, longwave radi-ation and precipitation; atmospheric CO 2 concentration; and solar radiation) during the period 1960-2006.Air temperature and longwave radiation were considered together, since downward longwave radiation can be expressed as a function of near-surface air temperature (e.g., Brutsaert, 1975).For air temperature and longwave radiation, we linearly regressed the annual values over time and removed cumulative changes due to the trend since 1960 from each subsequent year.For annual total precipitation and annual average shortwave radiation, we linearly regressed the annual values, computed each year's ratio of detrended to original annual values, and multiplied all original daily values by that ratio for each day within the year.For detrended atmospheric CO 2 , we used the 1960 concentration level for the entire period 1960-2006.Trends in the forcing variables were removed in cases when the trend was significant at the 0.05 level.At the 0.05 significance level, the entire domain experienced increasing trends in air temperature (0.0322 K yr −1 ), precipitation (0.5183 mm yr −1 ), [CO 2 ] (1.4009 ppm yr −1 ), and downward longwave radiation (0.0670 W m −2 yr −1 ), and a decreasing trend in downward shortwave radiation (−0.0385W m −2 yr −1 ), which is consistent with Fang and Yihui (2009;Table 2).

Sensitivities to climate drivers as a function of climate
We defined the sensitivity coefficients (α) of CH 4 emissions to long-term changes in the driver variables as the following partial derivatives:   Note: r 0 * is the reference CH 4 production rate per unit annual average LAI (r 0 * is related to the original r 0 parameter from Walter and Heimann (2000) by r 0 * = r 0 /LAI avg as described in Bohn et al., 2013); xvmax is the maximum CH 4 oxidation rate; rkm is the Michaelis-Menten constant for CH 4 oxidation; rq 10 is the Q 10 value for the CH 4 production rate; oxq 10 is the Q 10 value for the CH 4 oxidation rate; and tveg is a dimensionless integer value ranging from 0 to 15 that indicates the strength of plant-aided transport.
where the total change in annual methane emissions due to climate change CH 4 = α P × dP + α TLW × dT air + α CO 2 × dCO 2 + α SW × dSW + interaction.The CH 4 , T , [CO 2 ] and SW values in this relationship were annual average values, while P was annual total precipitation.We computed the sensitivity coefficients at each grid cell by first computing the time series of differences between the historical and control emissions and then performing a lin-ear regression between the differences in CH 4 and the differences between historical and detrended values of the driver variable.We then created maps of these sensitivities.To characterize the dependence of these sensitivities on climate, we divided the domain's grid cells into groups by their 46-year (1961-2006) average historical JJA T and JJA P , in increments of 2 • C and 20 mm, respectively (JJA T and P were chosen as independent variables for purposes of characterizing sensitivities rather than annual average T and P because the majority of annual CH 4 emissions occur during the growing season).Then, we computed the average sensitivities in each group and plotted them as a function of JJA T and P .This gave us two-dimensional matrices of sensitivities.Grid cells with the same JJA T and P conditions could come from quite different locations in the study domain; thus the resulting averaged sensitivities were not overly influenced by the characteristics of a single region.

Identifying the dominant emission controls
We calculated the correlation coefficients between the time series of CH 4 emissions and the various drivers at each grid cell, giving us a map of dominant controls (those with the highest correlations) across the domain.Similar to the sensitivities in Sect.2.4, we created two-dimensional matrices of correlations as a function of JJA T and JJA P .

Future projections
We generated two future projections of CH 4 emissions over the period 2007-2106: a process-based projection, in which we ran our modeling framework with future meteorological forcings, and a sensitivity-based projection, in which we applied the four sensitivity coefficients computed in Sect.2.4 to projected future forcings.To generate meteorological forcings for the future projections, we computed the monthly changes in meteorological forcings from the 4. CCSM4 (Community Climate System Model version 4) RCP4.5 (+4.5 W m −2 Representative Concentration Pathway) projection (which falls near the middle of the set of all CMIP5 (Coupled Model Intercomparison Project Phase 5) RCP4.5 projections) over the period 2007-2106, relative to the period 1996-2005, and applied these changes to the Sheffield et al. (2006) meteorology.
Based on the sensitivity matrices, and given a reference climate condition and corresponding CH 4 emission rate, we can derive the projected emission rate via CH 4 (t + 1) = CH 4 (t) + α P (T (t), P (t)) • (P (t + 1) where t is the year; CH 4 (t), T (t), P (t), and CO 2 (t) are the average values of annual CH 4 , JJA T , JJA P , and [CO 2 ] for the current grid cell over the last 10 years; and the coefficients are those defined in Eq. ( 1).Here we assume that interactions among the individual climate forcings are negligible.We check this assumption in Sect.3.2.1.We also assume that, as a grid cell's average T and P change, its sensitivities to drivers will evolve to resemble the current sensitivities of cells at the new (T , P ) coordinates.We discuss the validity of this assumption in Sect.4.2.

Parameter uncertainty
Our baseline simulations used the median parameter values of Bohn et al. (2013) as described in Sect.2.2.However, to assess the effects of parameter uncertainty on our results, we also generated an ensemble of 18 simulations using randomly sampled parameter values from the posterior distributions of Bohn et al. (2013; Table 1).The parameters that we examined included r 0 * (the reference CH 4 production rate per unit annual average LAI), xvmax (the maximum CH 4 oxidation rate), rkm (the Michaelis-Menten constant), rq 10 and oxq 10 (the Q 10 values for the temperature dependencies of the CH 4 production and oxidation rates, respectively), and tveg, a dimensionless integer value ranging from 0 to 15 that indicates the strength of plant-aided transport.The posterior distribution of tveg, which was held constant at a value of 12 in Bohn et al. (2013), was determined via Bayesian estimation from an ensemble of 3000 simulations that randomly sampled values of tveg across the range 0 to 15 and sampled values of all other parameters from their posterior distributions, in comparison with the observations of Glagolev et al. (2011).We did not vary the parameter pox, which represents the fraction of CH 4 oxidized in the root zone, as variations in tveg can compensate for variations in pox.Instead, we held pox constant at a value of 0.5, as in Walter and Heimann (2000) and Bohn et al. (2013).To account for uncertainty in our estimate of the border between the "south" and "north" regions, we performed two additional simulations, in which the entire domain used either the median "south" parameter set or the median "north" parameter set ("all-south" and "all-north", respectively).Adding these two simulations to our ensemble resulted in a total of 20 simulations.For each of these ensemble members, we constructed a distinct set of sensitivity matrices and created a sensitivity-based projection.

Historical simulation
Before examining simulated CH 4 emissions, we first evaluated model performance in simulating the environmental factors that are relevant to CH 4 emissions.The spatial distribution of simulated inundation extents was similar to that of the Surface Water Microwave Product Series ("SWAMPS") remote-sensing inundation product of Schroeder et al. (2010), with high concentrations in the WSL, Scandinavia, the HBL, and western Canada (Fig. 2a, b).VIC's inundated extent was biased low in western Canada, at about half the area given by SWAMPS.
To evaluate our simulated soil temperatures, we compared the distribution of continuous and discontinuous permafrost from the Circum-Arctic Map of Permafrost and Ground Ice Condition (CAMPGIC) map (Brown et al., 2014;Fig. 2c) with the VIC-simulated active layer depth (ALD) in the permafrost area (Fig. 2d).The spatial distribution of VIC's ALD was similar to the distribution of permafrost.An ALD of 1 m is an approximate threshold for "continuous permafrost" in the Brown et al. (2014) map.
We compared the simulated NPP distribution (Fig. 2e) with the MODIS MOD17A3 NPP product (Running et al., 2004;Fig. 2f).Model results and MODIS patterns matched reasonably well (spatial correlation 0.87), with a slight (about 6 %) overestimation of NPP in the boreal forest band between 55 and 65 • N latitude.
The spatial distribution of simulated average annual CH 4 emissions over the period 1960-2006 (Fig. 3) was similar to the distribution of wetlands (Fig. 1b), with notable concentrations in the WSL, Scandinavia, the HBL, and southern Canada.However, emissions were strongest in the boreal forest belt between 55 and 65 • N latitude, as a consequence of warmer temperatures, greater inputs of labile carbon (due to the higher rates of NPP there; see Fig. 2e, f), and the more productive "south" CH 4 parameter set that we used there.As an aside, the higher NPP values in the boreal forest belt do not necessarily imply that the peatlands there are forested, although some peatlands there do contain substantial numbers of trees (the VIC model does not distinguish between forested and non-forested peatlands).
We evaluated our simulated CH 4 emissions over three subdomains: the WSL, the HBL, and the high latitudes of the western hemisphere.Over the WSL, we compared our simulations with the estimate of Glagolev et al. (2011), which is based on in situ observations of mire landscape CH 4 emis- sions during 2007-2010 (Fig. 4).While our model tended to overestimate emissions in the middle of the domain, it captured the general north-south gradient in emissions.As to the total emission from the WSL area, Glagolev et al. (2011) estimated 3.91 ± 1.29 Tg CH 4 yr −1 , as compared with our estimate of 7.12 Tg CH 4 yr −1 .Our result here is also considerably higher than the estimate of Bohn et al. (2013) of 3.65 Tg CH 4 yr −1 , primarily because we (a) replaced that study's WSL-specific peatland maps (Sheng et al., 2004;Peregon et al., 2008) with the GLWD wetland map (Lehner and Döll, 2004), which attributes substantially higher wetland fractions to the region between 63 and 66 • N latitude than the WSL-specific maps do; (b) we replaced the WSL-specific assignment of "north" and "south" CH 4 parameter sets by the bioclimatic zone with the more general criterion of July LAI > 4 (Sect.2.2), which extended the region of more productive wetlands ("south" parameters) slightly further northward; and (c) used the meteorological forcings of Sheffield et al. (2006) instead of those of Adam et al. (2006).However, our estimate is within the range of estimates from inversions over the WSL, which range from 3.08 Tg CH 4 yr −1 (Kim et al., 2011)  Several studies have estimated total CH 4 emissions from all northern wetlands (Table 3), giving a range of 20-55 Tg CH 4 yr −1 over similar domains.Our model gives an estimate

Historical trends
Over the entire pan-Arctic domain, CH 4 emissions increased substantially over the period 1960-2006, with a trend of 0.158 Tg CH 4 yr −1 (Fig. 5a and Table 4, 4th column).Emissions from the control runs are shown in Fig. 5b-f.Defining the net impact of a driver as the difference between the historical trend in CH 4 emissions and the trend of the corresponding control run (Fig. 5g and Table 2, 4th column), we can see that air temperature and longwave radiation (TLW) had the largest impact on emissions (0.104 Tg CH 4 yr −1 , or 66 % of the historical trend), followed by CO 2 (0.030 Tg CH 4 yr −1 , or 19 %) and precipitation (0.015 Tg CH 4 yr −1 , or 10 %).The combined impact of TLW and P (TLWP), at 0.115 Tg CH 4 yr −1 , is slightly less than the sum of the impacts of TLW and P separately (0.119 Tg CH 4 yr −1 ), implying that these two drivers acted in opposition to each other to some extent but also indicating that the interaction between T and P was a relatively small effect.Locally, the effects of precipitation were often larger than those of CO 2 , but these effects largely canceled over the domain.

Sensitivity as a function of climate
The sensitivities of wetland CH 4 emissions to the climate factors we investigated varied in space or time and were strongly influenced by climate conditions.In Fig. 6a, which shows the distribution of spatial average annual CH 4 emissions as a function of 10-year average JJA T and P , "TLW" and "Tair LW" denote detrending of air temperature and associated downward longwave radiation; "CO 2 " denotes detrending of atmospheric CO 2 concentrations; "TLW+P" denotes detrending of both air temperature (and associated longwave radiation) and precipitation; "P " denotes detrending of precipitation; "SW" denotes detrending of downward shortwave radiation; and "inter" denotes the difference between "TLW" and "TLW+P".maximum CH 4 emissions occur along a "ridge" of slope 13 mm K −1 for JJA T values above 285 K and JJA P values above 120 mm.Consequently, increasing one factor (P or T ) while holding the other factor constant may cause CH 4 emissions to increase or decrease, depending on the current climate state of the wetland.Under relatively cold or dry conditions, emissions tend to increase with increasing T and P .However, at high P values, emissions decrease with increasing P , due to the inhibition of NPP under inundated conditions in the VIC model (Bohn et al., 2013).At high T values, emissions decrease with increasing T , due to increased oxidation of CH 4 as higher evaporation rates draw down the water table (Bohn et al., 2007).
Temporal correlations between historical annual CH 4 emissions and the three most important climate drivers (JJA T , JJA P , and JJA CO 2 ) were fairly consistent with this pattern (Fig. 6b-d).Correlations between annual CH 4 emissions and JJA T (Fig. 6b) were highest when JJA T is to the left of (colder than) the ridge of maximal emissions in Fig. 6a, and lowest (negative, in fact) to the right of (warmer than) the ridge.Similarly, correlations with JJA P were highest below (drier than) the Fig. 6a ridge and lowest (negative) above  (wetter than) the ridge, although this pattern broke down for JJA T below 285 K, where temperature limitation dominated the response and correlations with JJA P were only weakly positive or negative.Correlations with JJA CO 2 were moderately positive at all but the most extreme JJA T and P conditions, implying that CH 4 emissions generally benefit from CO 2 fertilization, via an increased input of carbon substrate into the soil.These differing responses of wetland CH 4 emissions to climate factors displayed strong geographic patterns, as a function of local climate (Fig. 7).In Fig. 7a, the ensemble median correlations between CH 4 and JJA T are represented on a blue (positive) to yellow (negative) color gradient.Similarly, correlations between CH 4 and JJA P are represented on a red (positive) to green (negative) color gradient.Therefore, blue indicates a strong positive temperature control on CH 4 emissions (T +), and this can be thought of as too cold for maximum emissions; yellow indicates a strong negative temperature (T −) control (too warm); green indicates a strong negative precipitation ("P −") control (too wet); and red indicates a strong positive precipitation ("P +") control (too dry).In general, northern cells were T + dominated (blue), due to the low summer air temperatures that they experience.These blue regions corresponded approximately to the distribution of permafrost (Fig. 1d).Moving southward, emissions became P + dominated (red).Southern West Siberia is relatively dry and warm, thus showing both P + and T + controls (orange).However, in the northernmost regions of Alaska and Canada (where inundation fractions were high, see Fig. 2b), we saw predominantly P − control (green).Comparison of this figure with Fig. 2b also shows that P + and T + (orange) areas were associated with smaller inundated area fractions and warmer temperatures, due to deeper water tables and greater oxidation rates.
Parameter-based uncertainties in the correlations (Fig. 7b), expressed as the range of correlations across the ensemble, were generally small (< 0.3) in both the T and P dimensions, except for P −-limited (green) regions in northeastern Canada and central Tibet and the northern portion of the T +limited region in north-central Canada.The general pattern of P + limitation in the southern reaches of the domain and T + limitation in much of the northern reaches of the domain appeared in all ensemble members.
Correlations between emissions and drivers tell us which driver is most influential at a given location.However, the sizes of the correlations are affected by both the relative sensitivities of emissions to the drivers and the relative amplitudes of the drivers' signals.It is therefore useful to consider the sensitivities alone.Sensitivities of annual emissions to the three main drivers (JJA T , JJA P , and JJA CO 2 ) were markedly higher outside the continuous permafrost zone than within it (Fig. 8).To first order, the explanation for this pattern is the general insensitivity of CH 4 emissions to all drivers at low temperatures, evident in Fig. 6a.Nevertheless, there were important differences among the distributions; for example, emissions in eastern Canada and eastern Siberia showed strong sensitivity to T but weak sensitivity to P and CO 2 .Spatial correlations between these sensitivities and various hydrologic and ecological terms, listed in Table 5, give some indication of which processes were most influential.The sensitivity of CH 4 emissions to JJA T (Fig. 8a) was most highly correlated (r = 0.30) with April-May snow water equivalent (AM SWE), which is consistent with a lack of water limitation, due to larger spring snowpacks leading to wetter summer conditions.Similarly, the sensitivity of emissions to P (Fig. 8b) was larger in absolute magnitude (positive or negative) where temperatures were warm, allowing for a higher (temperature-dependent) CH 4 production rate to be affected more dramatically by oxidation under drier conditions and reduced carbon input under wetter, more inundated conditions.The lack of strong correlations between the sensitivity to P and the various environmental factors in Table 5 may be the result of relatively high spatial heterogeneity in P and wetland moisture conditions (e.g., inundation), in comparison with those of T , leading to more "noise" in the relationships between them.Finally, the sensitivity of emissions to CO 2 (Fig. 8c) was most strongly correlated (r = 0.45) with NPP (Fig. 2f), which is consistent with the relationship between rates of carbon input into the soil and NPP in the model of Walter and Heimann (2000).Because relatively warm conditions and high NPP are associated with boreal forests, the geographic distributions of sensitivities to all factors also bore a strong similarity to the distribution of boreal forest.

Process-and sensitivity-based projections
To create a projection of future CH 4 emissions based on the climate sensitivities, we computed matrices of the sensitivity of aggregate emissions to each driver as a function of JJA T and P (Fig. 9a, c, e), similarly to the earlier correlation matrices (Fig. 6).To ensure that sensitivities exist for all possible future combinations of JJA T and P in the projection, we filled gaps in the matrices via a 3 row × 3 column window with a Gaussian kernel with σ = 1.Similar to the correlation matrices discussed in Sect.3.2.2, the sensitivities to JJA T , JJA P , and [CO 2 ] all exhibited maximum values along a diagonal "ridge" for T > 285 K and P > 120 mm (which correspond to the climate conditions in which boreal forest is found).For the sensitivities to JJA T and [CO 2 ], the ridges had similar slopes of approximately 30 mm K −1 .Sensitivities to JJA T were negative for P < 50 mm and 285 K < T < 291 K, due to increasing CH 4 oxidation above the water table with increasing temperature.In contrast, the ridge of maximum sensitivities to JJA P had a lower slope of about 12.5 mm K −1 , with a region of negative sensitivities for P > 190 mm and 287 K < T < 293 K, due to reduced productivity under inundated conditions.Again, sensitivities to all drivers were nearly 0 for JJA T < 285 K, due to the nonlinear temperature dependence of CH 4 production as well as the tendency for wetlands in that temperature range to be less productive (and therefore use the less productive "north" parameter set).Uncertainty in the methane model parameters (across the ensemble of random parameter combinations) led to a wide range of sensitivity values (Fig. 9b, d,  and f).However, the contours of the matrices of the individual ensemble members had similar shapes, so that regions of higher or lower sensitivities occurred in similar locations in climate space.The ensemble of sensitivity-based projections, created by applying these sensitivity matrices to meteorological forcings based on the CCSM4 RCP4.5 projection over the period 2006-2106, yielded a similar trajectory of CH 4 emissions to the projection from our process-based model (Fig. 10).Both the process-based (black) and median sensitivity-based (blue) projections agreed that emissions will initially remain relatively constant from 2007 to 2026 (in response to relatively little trend in air temperatures over the period; Fig. 10b) and then resume their increase.For the period 2056-2065, the process-and median sensitivity-based projections reached 46.1 and 43.4 Tg CH 4 yr −1 , respectively (132 and 124 %, respectively, of the 1997-2006 level).By the end of the century (2096-2105), they reached 50.1 and 48.3 Tg CH 4 yr −1 (142 and 138 %, respectively, of their 1997-2006 levels).Uncertainty in the methane model parameters led to a range of 39 to 57 Tg CH 4 yr −1 in sensitivity-based end-of-century emissions at the 95 % confidence level.However, the other members of the uncertainty ensemble followed trajectories that were similar to the median sensitivity-based projection over the course      While the two projections agreed on long-term behavior, their year-to-year variability disagreed at times, with the median sensitivity-based projection sometimes anticorrelated with the process-based projection.This is likely due to our construction of average sensitivities over all grid cells having similar climate conditions, which ignored the influence of local land cover, topography, and soils.Thus, during some years in some grid cells, our sensitivity matrices may have indicated a sensitivity of opposite sign to that of the processbased model, due to the grid cell's "ridge" of maximum emissions occurring in a different location in T −P space than in the domain-average matrix.Nevertheless, the general agreement in the long-term, domain-wide behavior implies that the sensitivity-based method captured the aggregate response of wetland CH 4 emissions to climate reasonably well.
Geographically, the regions of largest increases in emissions during the next century were in the boreal forest belt (Fig. 11a, c).This behavior was fairly consistent across the ensemble of methane parameter sets, with the exception of uncertainties > 30 % of the median in southern Canada and northwestern Siberia (Fig. 11b).These increases in emissions began at the southern edge of the domain and spread northward over time, corresponding to a northward shift in the types of controls exerted by climate factors, as shown in Fig. 12. Between 1997-2006and 2096-2105, areas of P + control (red and pink) migrated northward by 10-20 • of latitude, into territory that was previously under T + control (blue; Fig. 12, left).In other words, wetlands between 55 and 65 • N latitude that were previously colder than optimal experienced warming without a sufficient corresponding increase in precipitation, leading to their becoming drier than optimal and increasing their positive response to increases in precipitation.Other regions of historically T + control with large lake areas (e.g., Finland and northern Canada) were replaced by P − control (green) as they warmed.These patterns were robust across the parameter uncertainty ensemble, with large uncertainties primarily confined to northeastern Canada (Fig. 12 right).
To investigate the role of inundated area in the interannual variation of methane emissions, we calculated the changes in inundated fraction (Fig. 13a-c) and mean water table levels (Fig. 13d-f) of the process-based projection between the periods 1960-2006 and 2081-2100.The large areas of drying (reduced inundated area and falling water tables) in southern Canada and Alaska in Fig. 13c and f are consistent with the increase in the extent of P +-limited (red) wetlands in those same places over the 21st century, shown in Fig. 12.

Historical climate controls on CH 4 emissions
Our analysis indicates that summer air temperature increases explain almost two thirds of the long-term trend in CH 4 emissions over the last half century over the pan-Arctic domain.Precipitation had a smaller net effect (it explains only 10 % of the long-term trend), but this is due in part to spatial heterogeneity in the historical trends of P and their effects on CH 4 , leading to partial cancellation over the pan-Arctic domain.Nevertheless, the dominant role of air temperature in the pan-Arctic is not entirely surprising, given that the region is generally cold, leading to temperature limitation on metabolic rates.Our map of the historical controls on emissions (Fig. 7) corroborates this notion, since most of the region has historically been T + limited.This finding is largely consistent with Bloom et al. (2010), who also found that air temperature was the dominant factor controlling CH 4 emissions at high latitudes.However, our finding of strong P + limitation in the band between 50 and 60 • N (Fig. 7) is at odds with Bloom et al. (2010).This discrepancy may be due to a lack of variability in GRACE observations there (Bohn et al., 2015) or the inability of the global linear regression used by Bloom et al. (2010) to capture the location-and climate-dependent sensitivities accounted for by process-based models and the sensitivity-based approach that we have used here.
Within the pan-Arctic domain, we found strong geographic patterns in climate controls on CH 4 emissions.Sim-ilar (observation-rather than model-based) analyses have been performed on NPP (Nemani et al., 2003) and ET (Teuling et al., 2009).Our study shares some similarity in conclusions.For example, these studies show that CH 4 , NPP and ET are all T + controlled around Hudson Bay and in Scandinavia, and P + controlled in the wetlands of southwestern Canada.This is not surprising, because NPP and ET are both tightly linked with CH 4 production: NPP determines how much carbon can be converted to CH 4 , while ET is positively correlated with soil moisture content, as is the CH 4 emission rate.In the WSL, the wetlands in the south are P + and T − controlled, suggesting that this area is much drier than the north, with more CH 4 emitted as the water tables are drawn down during summer (Bohn et al., 2007).NPP in this southern area is in transition from T limited to P limited (Nemani et al., 2003), which is consistent with CH 4 .In a recent process-based study, Liu et al. (2015) also found that ET in southern Siberia is P + limited.
Despite their similarities, there are some differences in the spatial distributions of controls between our and previous studies.In Nemani et al. (2003), NPP over northern Europe and West Siberia is almost entirely limited by temperature and radiation, while in our results, CH 4 is P + limited over a considerable area.This is due in part to the nearly negligible role shortwave radiation plays in CH 4 emissions (Fig. 5), in part to the drier optimal soil moisture conditions for upland vegetation (included in the Nemani et al.NPP analysis), relative to wetland plants (which we focus on here), and in part to the rapid drop in CH 4 emissions as the water table is drawn down beyond a few centimeters.Similarly, the area of P + limitation of ET in western Canada in Teuling et al. (2009) is smaller than the area of P + limitation of CH 4 emissions in our study.This can also be explained by the presence of forested uplands in this area, where the moisture deficit in upper soil layers from low precipitation is partly compensated for by water extracted from deeper soils.Thus, only those places with a considerable shortage of water will show up as P + in the Teuling et al.ET map.
The validity of our results depends on our model's temporal behavior, which is subject to both model uncertainty and parameter uncertainty.Verification opportunities include in situ observations and atmospheric model inversions.Both are problematic, due to the paucity of long in situ observational records in the first case and the errors to which inversions are subject in the second (see Bohn et al., 2015).Nevertheless, in Bohn et al. (2015), the interannual variability of our modeling framework (called "UW-VIC" therein) was assessed over the West Siberian Lowland over the period 1993-2004, relative to observations, several atmospheric inversions (including those of Bousquet et al., 2011), and many other processbased models.While there was little agreement across these data sets in terms of interannual variability, the process-based models (including UW-VIC) that employed soil physics formulations appropriate to high latitudes and realistic relationships between CH 4 emissions and water table depth tended to   be more similar to the inversions than those that did not.Our investigation of parameter uncertainty (Figs. 9 and 10) revealed a substantial range in sensitivities and end-of-century CH 4 emissions but made little difference to the shape of the trajectory over the next century or the spatial distribution of climate controls.Thus, we believe our findings here are robust with respect to parameter uncertainty.However, investigation of the impacts of model uncertainty on climate controls on CH 4 emissions using other model formulations would be useful.

Sensitivity-based future projections
Our sensitivity estimates provide a simplified description of wetland behavior and is in effect, a linearization of our process-based model.Nevertheless, the similarity between our process-based and sensitivity-based projections suggests that our domain-averaged sensitivities capture most of the dependence of CH 4 emissions on climate conditions, as represented within our modeling framework.Our projected emissions are comparable to those of other process-based studies.Our estimate of a 24-32 % increase in pan-Arctic CH 4 emissions by mid-century is comparable to the 25 % increase estimated by Anisimov (2007).Over northern Eurasia, our estimate of end-of-century emissions is 21.5 Tg CH 4 yr −1 , similar to the estimate of 25.1 ± 3.7 Tg CH 4 yr −1 by Zhu et al. (2011).The widespread warming and drying of wetlands and consequent reduced sensitivity of emissions to warming in our projections are consistent with similar findings in other studies (Koven et al., 2011;Riley et al., 2011;Ringeval et al., 2011;Lawrence et al., 2015).
Our characterization of the sensitivities of emissions to climate requires the assumption that, as a grid cell's climate changes, its future sensitivities will come to resemble those of cells with similar climate today, in essence attributing climate sensitivities completely to current climate state.Several studies have, however, found associations between vegetation and CH 4 emissions (Glagolev et al., 2011;Lupascu et al., 2012;Levy et al., 2012;Olefeldt et al., 2013).In particular, Olefeldt et al. (2013) found that emission rates from sedge-dominated wetlands are not only higher but also more sensitive to changes in both soil temperature and water table depth than are emission rates from non-sedge-dominated wetlands.On the other hand, dynamic vegetation models suggest that vegetation communities will migrate northward with future climate change (e.g., Kaplan and New, 2006;Alo and Wang, 2008), potentially bringing with them any characteristics (e.g., aerenchyma) that enhance CH 4 emissions.To the extent that vegetation communities can migrate in step with climate change, our sensitivity matrices would still be applicable.Nonetheless, this suggests an interesting avenue for future research.

Future changes in the dominant controls
In our future projections, we found that much of the region will shift from T + limitation (colder than optimal) to T − and P + limitation (warmer and drier than optimal).This large-scale shift towards the warm and dry side of the "ridge" www.biogeosciences.net/12/6259/2015/Biogeosciences, 12, 6259-6277, 2015 of maximum emissions implies that air temperature will play a smaller role in end-of-century emissions than at present, for two reasons: first, the positive response to an increase in temperature in the northern portion of the domain will be partially or completely canceled by the negative response from the southern portion, and, second, the response to precipitation will increase due to the widespread drier-than-optimal conditions.This suggests that, beyond the year 2100, emissions may level off or even decrease under further climate change, unless precipitation can increase sufficiently to compensate for the increases in air temperature.The larger future role of P in controlling pan-Arctic CH 4 emissions may lead to greater uncertainty in future projections beyond 2100, due to the poorer performance and greater lack of agreement of global climate models in projecting future precipitation than temperature (Hawkins and Sutton, 2011;IPCC, 2013).
There are additional reasons to think that T will play a reduced role in the future.There is some indication that the metabolic impacts of higher temperatures have been overestimated by most models, as most studies neglect acclimatization.Koven et al. (2011), for instance, found that soil microbial communities essentially adapt to warmer soil temperatures and CH 4 emissions rates return to their previous levels.Koven et al. (2011) showed that acclimatization could eliminate over 50 % of the increase in emissions over the pan-Arctic by the end of the century that would otherwise occur.Under such conditions, the primary effects of increased T would then be on drying out the wetlands through increased ET.In addition, because our model did not simulate dynamic vegetation phenology, we did not account for increased transpiration arising from CO 2 fertilization, which also would have a drying effect on the wetlands (the wetland-climate CH 4 feedback as discussed by Ringeval et al., 2011;Koven et al., 2011;and Stocker et al., 2013).Including these effects in the model on which our sensitivities were based would likely reduce the sensitivity of future emissions to further increases in T and perhaps even change the sign of the sensitivity to negative in some water-limited locations.Thus, our estimates of the expansion of the water-limited zone and the reduction of the role of T may be considered a lower bound.

Conclusions
We performed an historical simulation of wetland CH 4 emissions for the pan-Arctic domain for 1948-2006.In addition, we performed five experiments that investigated the sensitivities of CH 4 emissions to changing climate and two future projections over the period 2007-2106 -one process-based and the other based on CH 4 emission sensitivities to T , P and CO 2 .Our main conclusions are as follows: 1. We estimate the annual CH 4 emissions from pan-Arctic wetlands averaged over 1997-2006 at 35.0 ± 6.7 Tg CH 4 yr −1 .This is within the range of previous estimates but somewhat toward the higher end.
2. Based on our model, climate change over the last half century has led to a substantial (20 %) increase in total emitted CH 4 , with increases in air temperature (and associated downward longwave radiation) being the dominant driver.Increases in temperature and [CO 2 ] were responsible for over 84 % of the inferred increase in emissions.Most of the remainder is attributable to changes in shortwave radiation (decreasing) and precipitation (increasing).
3. The dominance of air temperature is corroborated by the predominance of temperature-limited wetlands throughout most of the domain, with water-limited wetlands primarily occupying only the southernmost portion of the domain (south of 60 • N latitude).
4. Both process-based and sensitivity-based projections agreed that wetland CH 4 emissions from pan-Arctic wetlands will increase to 138-153 % of present-day levels by the end of this century.Because this study did not account for potential acclimatization or the wetlandclimate CH 4 feedback resulting from CO 2 fertilization, this projected increase may be overestimated.
5. As future climate across the pan-Arctic becomes warmer, northern wetlands are likely to shift from the current temperature-dominated state toward a more precipitation-dominated state due to a lack of sufficient increase in precipitation to compensate for higher evapotranspiration and resultant soil drying.The resulting sensitivity of CH 4 emissions to further warming may then level off or even become negative.
Author contributions.All coauthors jointly conceived and designed this study.X. Chen performed all model simulations and data analysis.X. Chen prepared the manuscript with contributions from all coauthors.

Figure 3 .
Figure 3. Average annual CH 4 emissions over the study domain for 1960-2006.

Figure 5 .
Figure 5.Time series of domain-averaged annual methane fluxes from (a) the historical simulation; (b-f) the five climate control runs, in each of which one climate driver was detrended starting in 1960; (g) differences between historical simulation in (a) and the control runs (b-f)."TLW"and "Tair LW" denote detrending of air temperature and associated downward longwave radiation; "CO 2 " denotes detrending of atmospheric CO 2 concentrations; "TLW+P" denotes detrending of both air temperature (and associated longwave radiation) and precipitation; "P " denotes detrending of precipitation; "SW" denotes detrending of downward shortwave radiation; and "inter" denotes the difference between "TLW" and "TLW+P".

Figure 6 .
Figure 6.Panel (a): the 1960-2006 average annual CH 4 emission over JJA (June-July-August) T and JJA P space; panels (b-d): correlation between 1960-2006 annual CH 4 emission and JJA drivers in the same T −P space.

Figure 7 .
Figure 7. Spatial distributions of ensemble median (left) and range at 95 % confidence level (right) of correlations between annual CH 4 emissions and JJA T and P .The green-red and yellow-blue axes depict the strength of correlation (−1 to 1) with JJA P and JJA T , respectively.

Figure 8 .
Figure 8. Spatial distributions of sensitivities of CH 4 to climate drivers.Panel (a): sensitivity to air temperature; panel (b): sensitivity to precipitation; panel (c): sensitivity to [CO 2 ].

a
JJA T : June-July-August average air temperature; b JJA P : June-July-August total precipitation; c JJA F inund : June-July-August inundated area fraction; d AM SWE: April-May average snow water equivalent; e JJA LAI: June-July-August average leaf area index; f ALD: maximum annual active layer depth; g Annual NPP: annual net primary productivity; h JJA Tsoil: June-July-August average temperature in the top 10 cm of the soil column.i Extreme values of sensitivity (> 0.005 g CH 4 m −2 yr −1 per mm change in JJA P ) were ignored; these occurred in 164 cells, out of 2049 cells containing wetlands.

Figure 9 .
Figure 9.The 1960-2006 average T , P and CO 2 sensitivities (a, c, and e, respectively) of CH 4 emissions in JJA T and JJA P space using the median methane model parameter set and their ranges at the 95 % confidence level across randomly sampled methane model parameter sets (b, d, and f, respectively).

Figure 10 .
Figure 10.Historical and projected annual methane emissions and climate drivers over the pan-Arctic from 2007 to 2106.Panel (a): sensitivity-and process-based projections (blue and black, respectively) of methane emissions from northern wetlands during 2007-2106, with historical simulation (red) 1948-2006.Parameter-based uncertainties in the sensitivity-based projection are plotted as the yellow and green envelopes (50 and 95 % confidence bounds, respectively); panels (b-d): climate conditions for projections.The end-of-century window for time slice analysis (2096-2105) is denoted with vertical solid and dashed lines in (a).

Figure 11 .
Figure 11.Ensemble median (a) and range (b) of average annual end-of-century (2096-2105) CH 4 emissions for the sensitivity-based projection and the difference between the median and the annual emissions of year 2006 (c).

Figure 12 .
Figure 12.Spatial distributions of ensemble median (left) and range at 95 % confidence level (right) of correlations between annual CH 4 emissions and JJA T and P for the period 2081-2100 of the sensitivity-based projection.The green-red and yellow-blue axes depict the strength of correlation (−1 to 1) with JJA P and JJA T , respectively.

Figure 13 .
Figure 13.Changes in inundated area fraction and water table position during the historical period (1960-2006) and projection period (2081-2100) for the process-based projection.Panel (a) is the average inundation fraction during 1960-2006, (b) is the average of 2081-2100, (c) is the difference between these two averages (b-a).Panels (d-f) are similar for water table positions.

Table 2 .
Trends in spatial average climate factors from 1960 to 2006.

Table 3 .
Estimates of total CH 4 emissions over the study domain.Variable Infiltration Capacity plus Terrestrial Ecosystem Model, b Model for Atmospheric Transport and Chemistry, c Wetland Methane Emission Model. a

Table 4 .
Trends in CH 4 emissions from historical and control simulations from 1960 to 2006.All values are in units of Tg CH 4 yr −1 .

Table 5 .
Spatial correlation coefficients between sensitivities and environmental factors.