The effects of surface moisture heterogeneity on wetland carbon fluxes

Introduction Conclusions References


Introduction
As the world's largest natural source of methane (CH 4 ), wetlands are an important component of the global carbon cycle (Fung et al., 1991).Northern wetlands additionally have been a large net sink of carbon dioxide (CO 2 ) since the last Ice Age (Smith et al., 2004).Because the fluxes of both of these greenhouse gases are highly sensitive to soil moisture and temperature, there is concern that wetland CH 4 and CO 2 emissions could provide a positive feedback to climate change (Ringeval et al., 2011;Wania et al., 2009;Tang et al., 2008).Of particular concern are carbon-rich boreal and arctic wetlands, which comprise up to 50 % of the total global wetland area (Lehner and D öll, 2004) and occupy one of the world's hotspots of historic and projected climate change (Serreze et al., 2000;Diffenbaugh and Giorgi, 2012).Despite the importance of these ecosystems to the global carbon cycle, substantial uncertainties remain in estimates Introduction

Conclusions References
Tables Figures

Back Close
Full of their current and future greenhouse gas emissions, due to uncertainties in both their emissions per unit area and the extents of their contributing areas.Soil moisture's control on wetland carbon fluxes is complex and spatially heterogeneous.Because anoxic conditions in water-saturated soil promote methanogenesis and inhibit aerobic respiration, the water table position controls the partitioning of decomposition into CO 2 and CH 4 (Walter and Heimann, 2000) and, therefore, decomposition's net greenhouse warming contribution.In northern peatlands, photosynthesis is also inhibited when the water table reaches the surface (Frolking et al., 2002;Wania et al., 2009).Within boreal wetlands, a mosaic of lakes (Repo et al., 2007;Smith et al., 2005), seasonally-varying inundation (Schroeder et al., 2010), and spatially-varying water table depth on the scale of meters (Eppinga et al., 2008) gives rise to a saturated zone, whose area varies in time, and where CH 4 fluxes are at their maximum (Bohn et al., 2007).
Because the dependence of carbon fluxes on soil moisture is nonlinear, it is important to account for their spatial distribution in estimating these fluxes (Baird et al., 2009), especially their response to climate change (Bohn and Lettenmaier, 2010).To date, large-scale models have not accounted for all of the effects of heterogeneity on carbon fluxes, nor have they agreed on consistent definitions of the fluxes' contributing areas (Melton et al., 2013).As described in Bohn and Lettenmaier (2010), most models have employed one of two schemes to describe wetland surface and subsurface water: either a uniform water table scheme (which ignores the saturated zone as well as the distribution of water table depths within the wetland) (e.g.Zhuang et al., 2004) or a wet-dry scheme (which neglects the contribution of the unsaturated wetlands to methane fluxes) (e.g.Ringeval et al., 2010;Gedney et al., 2004).When applied to a typical boreal wetland, these schemes' predictions of the response of methane emissions to future climate change can be subject to biases of up to ±30 %, compared to a distributed water table scheme (Bohn and Lettenmaier, 2010).At large scales, different assumptions about local methane emissions and contributing areas can result in large differences in the spatial distribution of methane emissions (Petrescu et al., 2010;Introduction Conclusions References Tables Figures

Back Close
Full  Melton et al., 2013).However, the impact of moisture heterogeneity on the total carbon budget at regional and global scales is still poorly known.Given the sensitivity of wetland carbon fluxes to sub-grid heterogeneity in surface and sub-surface water, our goal was to assess the influence of this heterogeneity on the carbon budget of boreal wetlands over the last half-century.We focused our study on the West Siberian Lowland, where recent intensive field campaigns across large areas (Sheng et al., 2004;Peregon et al., 2009;Glagolev et al., 2011) enabled us to evaluate the performance of models in reproducing the large-scale spatial patterns of carbon fluxes.We used a combination of large-scale models, remote sensing, and in situ measurements to examine the spatio-temporal distributions and carbon fluxes of the region's saturated and unsaturated wetlands.

Study domain
The Western Siberian Lowland (WSL, Fig. 1) is a low-lying region bounded by the Ural Mountains to the west, the Yenisei River and Central Siberian Plateau to the east, the Arctic Ocean to the north, and the Central Asian steppe to the south (Kremenetski et al., 2003).Containing between 590 000 and 680 000 km 2 of wetlands, and 70 Pg of carbon in its peat soils (10-15 % of the world's boreal carbon reservoir), the WSL is the world's largest high-latitude wetland region (Kremenetski et al., 2003;Sheng et al., 2004;Peregon et al., 2009).As shown in Fig. 1, permafrost occurs in the northern half of the domain (north of about 61 • N).Interspersed among the region's wetlands are thousands of lakes and ponds, ranging in morphology from bog pools to thaw lakes (Eppinga et al., 2008;Repo et al., 2007;Smith et al., 2005;Lehner and D öll, 2004).
Portions of the wetlands of the WSL undergo seasonal inundation after the region's snow melts in late spring and early summer (Schroeder et al., 2010).Introduction

Conclusions References
Tables Figures

Back Close
Full

Modelling framework
To account for all of the major components of the wetlands of the WSL accurately (lakes, inundated wetlands, and exposed wetlands, as depicted in Fig. 2), we employed a modified version of the Variable Infiltration Capacity (VIC) model (Liang et al., 1994), version 4.1.2as the land surface component of the modeling framework.The relevant features of VIC 4.1.2include simulation of permafrost using the scheme of Cherkauer and Lettenmaier (1999), thermal properties of organic soils (taken from Farouki, 1981), and a dynamic lake/wetland model that accounts for impoundment of surface water, thermal stratification, mixing, and ice cover (Bowling and Lettenmaier, 2010).Special enhancements include computation of carbon cycle processes such as net primary productivity (NPP) and soil respiration (Rh), with inhibition of both of these fluxes under saturated conditions; a distributed water table scheme; and a parameterization of wetland microtopography.The modified VIC model was linked to the wetland methane emissions model of Walter and Heimann (2000) for calculation of methane fluxes, as described in Bohn et al. (2007).In addition, we made a small modification to the Walter and Heimann (2000) methane emissions model to allow sensitivity to spatial variation in NPP.These features are described in detail in Appendices A1 and A2.

Meteorological forcings and spinup
Daily meteorological forcings for the period 1948-2010 were created by rescaling the NCAR/NCEP reanalysis-derived fields of Sheffield et al. (2006) to match the monthly statistics of gridded monthly observations (precipitation and wind speed from Wilmott and Matsuura, 2001; air temperature from Mitchell and Jones, 2005), with corrections to precipitation for gauge undercatch (Adam and Lettenmaier, 2003) and orographic effects (Adam et al., 2006).Daily short-and longwave radiation and humidity, and hourly values of all meteorological forcings, were derived using methods described in Bohn et al. (2013).Introduction

Conclusions References
Tables Figures

Back Close
Full To achieve realistic initial soil temperatures without a lengthy spin-up, soil temperatures were initialized with soil temperature profiles from Troy et al. (2011).This initialization was followed by a second 30 yr spinup using forcings from the period 1948-1957 three times.Soil carbon pools required a different approach, because they are likely only 50-90 % of their long-term equilibrium size, assuming an average wetland age of 10 000 yr (Sheng et al., 2004) and drawing from Fig. 6 of Wania et al. (2009).Therefore, we first iteratively ran our soil respiration model offline over the period 1947-1958 to find the long-term equilibrium carbon pool densities (defined by a net change in carbon pool storage of less than 0.1 g C m −2 over the 10 yr period).Then we rescaled the carbon densities by a constant factor (equal to 0.818 for simulations using optimal parameters) across the entire domain, so that the total storage in 1948 equaled the estimate of 70.2 Pg C from Sheng et al. (2004).

Hydrology
As shown in Fig. 3, we divided the domain into 278 cells using the EASE-Grid 100 km polar azimuthal equal area grid (Brodzik and Knowles, 2002).To account for different soil textures and flow regimes in lake-wetland systems versus uplands, we further divided these cells into dry and wet portions.The wet portion of each cell consisted of the cell's peatlands, taken from the database of Sheng et al. (2004); permanent lakes, taken from the Global Lake and Wetland Database (GLWD; Lehner and D öll, 2004); and wet tundra land cover classes, taken from Bartalev et al. (2003).Each cell's wet portion was modeled as a single lake-wetland tile, underlain by peat soil, with peat depth taken from Sheng et al. (2004) and hydraulic properties taken from Letts et al. (2000).Exposed land within the lake-wetland tile was assigned to one of three VIC land cover classes (bog, forested bog, or evergreen needleleaf forest), according to the proportions of forest and non-forest pixels from Bartalev et al. (2003).Parameters for these classes are listed in Table 1.LAI values were prescribed by the monthly Introduction

Conclusions References
Tables Figures

Back Close
Full average over all pixels from the MODIS LAI product (Myneni et al., 2002) within the wet portion of each cell.Because this study focused on wetland behaviors, all descriptions hereafter refer only to the wet portions of the domain.Because VIC's NPP and Rh are inhibited under saturated conditions and reduced to zero under inundated conditions (when low-lying vegetation is completely submerged), it was important to constrain the fractional areas of inundation and saturation (A inund and A sat , respectively).To do so, we used two remote sensing products: the coarse-resolution (25 km) AMSR-E/QuikSCAT-based passive-active microwave global inundation product of Schroeder et al. (2010), extended to all snow-free days from July 2002 through 2009 (available at http://wetlands.jpl.nasa.gov);and high-resolution (30 m) ALOS/PALSAR active microwave imagery from four main focus regions (delineated by boxes in Figs.3b and d) on several dates from 2006 and 2007.The PALSAR pixels were classified into various categories, including open water (no emergent vegetation) and wet soil (emergent vegetation, with the water surface ranging from a few cm above the soil surface to just below the soil surface), via the random forest method described in Whitcomb et al. (2009).The extent of the saturated zone was computed as the sum of open water and wet soil extents.Comparison of the PALSAR classifications to the Schroeder et al. product over each 100 km ease grid cell where they overlapped (Fig. 3a) indicated a good match between the Schroeder et al. product and the PALSAR open water extents (bias 0.0172, R 2 of 0.79).Therefore, we chose to treat the Schroeder et al. product as a measure of A inund .
For VIC's dynamic simulation of A inund , we created a two-part depth-area relationship for each grid cell.The first part, which describes the composite bathymetry of permanent water bodies, was computed by assigning a depth to each of the cell's lakes, as a function of the lake's size, and computing the cumulative surface area of lakes deeper than each of several pre-selected depths.Lake depths were computed via a linear regression of log(depth) vs. log(area) using data reported for bog pools (Belyea and Lancaster, 2002;McEnroe et al., 2009), Alaskan thaw lakes (Bowling et al., 2003), and larger West Siberian lakes from the International Lake Environmental Committee Introduction

Conclusions References
Tables Figures

Back Close
Full World Lake Database (ILEC, 1988(ILEC, -1993)).The areas of this depth-area relationship were then rescaled so that the total lake area fraction A lake matched the minimum monthly average value of A inund given by the Schroeder et al. product.The second part of the depth-area relationship described the additional surface area (A inund between A lake and A wet ) covered by seasonal inundation, as a function of excess water depth above the lake surface.This part of the relationship was generated by integrating the distribution of surface slopes from either the Shuttle Radar Topography Mission (SRTM) DEM (Farr et al., 2007) south of 60 • N , or the Advanced Spaceborne Thermal Emission and Reflection (ASTER) DEM (NASA, 2001) north of 60 • N .
The remaining lake parameter wfrac (outlet width as a fraction of the lake perimeter) was calibrated separately at each grid cell to minimize the bias between local simulated and observed JJA average A inund over the period 2002-2007.JJA average simulated A inund over the period 2001-2010 is shown in Fig. 3b.Within the wetland, two primary parameters needed to be determined: the fraction of the wetland covered by ridges (f ridge ) and the maximum subsurface flow rate (Dsmax).We held these parameters constant across the domain.We sampled f ridge and Dsmax uniformly across the ranges of their plausible values: 0.4 to 0.6 for f ridge , based on the findings of Peregon et al. (2009); and 0.0 to 0.8 mm day −1 for Dsmax, with the upper bound determined empirically by the maximum value of Dsmax for which A sat > A inund .We designated as "optimal" the values that approximately minimized the total bias between simulated and observed (via PALSAR) A sat across the locations and times of PALSAR observations.Optimal values for both f ridge and Dsmax equaled 0.

Biogeochemistry
Where possible, parameters for photosynthesis were taken from the BETHY model (Knorr, 2000).Values of the parameters V m and J m chosen for the new bog and forested bog land cover classes are listed in Table 1.Most parameters for aerobic soil respiration and methane emissions were taken from Sitch et al. (2003) and Walter and Heimann (2000), respectively.Three parameters associated with photosynthesis and soil respiration required calibration: the wetland photosynthesis inhibition parameter f inhib (Eq.A4), the soil respiration inhibition factor r sat (Eq.A6), and the soil respiration scaling factor k (Eq.A5).These parameters were held constant across the domain.We calibrated these parameters using CO 2 flux observations for the year 2000 from the Zotino Bog flux tower (Arneth et al., 2002), denoted with a yellow star in Fig. 1.We assessed the joint likelihoods of 1000 uniformly-sampled parameter combinations with respect to two objective functions: (1) the joint likelihoods of the time series of simulated vs. observed 10 day mean net ecosystem exchange (NEE); and (2) the mean soil carbon density of 90 ± 30 kg m −2 for the region around the flux tower, based on Sheng et al. (2004).Percentiles of the posterior distributions of parameter values are listed in Table 2.For the Walter and Heimann wetland methane emissions model, we calibrated five parameters: r0, xvmax, rkm, rq10, and oxq10.We used in situ observations of methane flux, soil temperatures, and water table depth from over 750 locations across West Siberia for the period 2006-2010 (Glagolev et al., 2011(Glagolev et al., , 2012;;Sabrekov et al., 2012), most of which were only monitored for a period of 1-2 days.To generate simulated fluxes at each site, soil temperatures, NPP, and water table distributions simulated by VIC for the grid cell containing the site were used as inputs to the methane emissions model.To form an objective function, we grouped the observations into five regional groups (depicted in Fig. 1).Within each regional group, we divided the observations into four water table depth categories (shown in Fig. 5) and computed the mean of the log-transformed simulated and observed fluxes for each category.For each regional group, our objective Introduction

Conclusions References
Tables Figures

Back Close
Full function was the joint likelihood of the categories' means.We found that using a single parameter set for all observation groups resulted in substantial positive biases in the northern three groups (groups 1-3 in Figs. 1 and 5).Therefore, we calibrated groups 1-3 and 4-5 separately, resulting in separate north and south parameter sets.Posterior distributions of values for these two parameter sets are listed in Table 2.We applied the south parameter set to grid cells having July wetland LAI values ≥ 2.5, and the north parameter set to all other cells.The boundary between these two sets of cells fell approximately at 61 • N, roughly coinciding with both the northern boundary of the middle Taiga (Glagolev et al., 2011) and the southern limit of permafrost.The resulting distributions of observed and simulated soil temperatures and CH 4 fluxes are shown in Fig. 5. Simulated soil temperature profiles (left column) matched observations reasonably well across all groups, except for a warm bias of 1-2 • C from 0 to 20 cm depth in the north (groups 1-2).Simulated CH 4 fluxes (right column) matched observed fluxes reasonably well in the south (groups 4 and 5), but less so in the north.In particular, the simulations substantially overestimated fluxes under saturated conditions at sites in the central WSL (groups 3-4), and substantially underestimated fluxes in the far north (group 1).

Historical simulations
Historical simulations covered the period 1948-2010, using the optimal (median posterior) values for all calibration parameters.To assess parameter uncertainty, we then randomly sampled calibration parameter values to generate 1200 simulations at a randomly-selected subset of 50 cells across the domain, the results of which were spatially interpolated to the other cells.To assess the influence of saturated area on carbon fluxes, we performed an additional control simulation employing a uniform water table scheme, in which each grid cell's water table depth distribution was replaced by its spatial average and fractional inundated and saturated areas were set to 0.

Present-day extents and carbon fluxes
In terms of both area and total carbon fluxes, unsaturated wetlands were the dominant hydrologic zone in the WSL.Estimates of the current JJA average extents (chosen because the vast majority of carbon fluxes occur during the summer) of the wetland zones over the period 2001-2010 are listed in Table 3.The total JJA saturated area (excluding permanent lakes) had a median value of 235 000 km 2 (with 1st and 99th percentiles of 165 000 and 332 000 km 2 , respectively).This represents 34 % (24 to 47 %) of the non-lake wetland area of the WSL (and 13 % of the total area of the WSL).
Inundated wetlands occupied 25 250 km 2 , or 4 % of the non-lake wetland area.
The unsaturated zone dominated the area of the WSL wetlands, as well as the region's carbon budget.As shown in Table 4, the total annual wetland NPP and Rh across the WSL are 163 and 149 Tg C yr −1 , respectively, and the contribution from unsaturated wetlands (129 and 119 Tg C yr −1 , respectively) comprises about 80 % of the total flux in both cases.Similarly, unsaturated wetlands emit 64 % of the total methane flux of 3.6 Tg CH 4 yr −1 .Because NPP exceeds Rh, the WSL wetlands as a whole are a net carbon sink of 11.5 Tg C yr −1 (sinks are listed as negative in Table 4); unsaturated wetlands account for 8.0 Tg C yr −1 , or 70 % of the region's wetland carbon sink.Assuming that CH 4 has a greenhouse warming equivalence of 21 over the span of a century (IPCC, 2007), CH 4 emissions dominate the region's total greenhouse warming potential (GHWP), leading to a total GHWP of 24.7 Tg CO 2 yr −1 ; unsaturated wetlands account for 56 % of the total.The dominance of unsaturated wetlands in the region's NPP and Rh may not come as a surprise, given that saturated wetlands should exhibit lower rates for these fluxes (rates per unit area in saturated wetlands tended to be 50-60 % of their rates in adjacent unsaturated wetlands).In contrast, the dominance of unsaturated wetlands in CH 4 BGD 10,2013 The effects of surface moisture heterogeneity on wetland carbon fluxes Introduction

Conclusions References
Tables Figures

Back Close
Full and GHWP fluxes is counterintuitive, because CH 4 fluxes per unit area tended to be much larger (by about 40 %) in the saturated zone.A second surprising result was that the control simulations, for which non-lake saturated area was always 0, yielded similar results to the primary simulations: an increase of only 5 % in NPP, Rh, and C net , and a decrease of only 2 % in CH 4 and GHWP, over the totals from the distributed scheme.
The large area of the unsaturated zone was partly responsible for these results.

Spatial distributions
Not only did unsaturated wetlands occupy a larger total area than saturated wetlands, but also the majority of the WSL's carbon fluxes occurred in the drier parts of the domain, where the saturated area fraction is small.To illustrate this, maps of the simulated 2001-2010 annual average states and fluxes (expressed as fluxes or densities per unit area of non-lake wetland) are shown in Fig. 5. JJA average A sat (Fig. 5a) is concentrated in the central and northern portions of the domain (north of 60 • N).
In contrast, the spatial distributions of NPP, Rh, and soil carbon density (C dens ), shown in Figs.5b-d, were all anti-correlated with A sat (correlation coefficients of −0.72, −0.71, and −0.64, respectively).Control runs using a uniform water table scheme (not plotted) yielded similar spatial distributions to the distributed water table scheme.Thus, the spatial distribution of these terms appears to have been primarily a reflection of the spatial distribution of VIC's MODIS-prescribed LAI, which had a correlation coefficient with A sat of −0.71 (although this distribution of LAI may itself be the cumulative result of thousands of years of NPP inhibition).
The spatial distributions of CH 4 , C net , and GHWP (Fig. 5e-g) demonstrated weaker dependence on A sat (correlations of −0.43, +0.30, and −0.30, respectively) but even in these cases, fluxes above the median rate tended to fall outside the region of greatest A sat .In particular, the CH 4 emissions (Fig. 5e) were substantially higher south of 61 between 57 • and 66 • N. Still, the majority of uptake rates greater than the median rate of 20 g C m −2 yr −1 fell in cells where A sat was less than 0.5.The distribution of GHWP (Fig. 5 g) was more complex, with the large C net uptake rates in the West driving GHWP negative, and the large CH 4 emissions rates in the South and East driving GHWP positive.Still, the largest positive or negative fluxes once again occurred outside the wetter center of the WSL.
Our simulated C dens distribution was somewhat corroborated by the distribution of observed soil carbon densities reported by Sheng et al. (2004) (Fig. 6).With the exception of a mismatch in the North (66-70 • N, 75-90 of the control simulations (Fig. 7d-f) led to a substantially greater proportion of emissions in the northern half of the domain.Thus, to avoid large CH 4 emissions in the central and northern WSL, it appears that simulations must account for both the lower methanogenesis rates north of 61 • N and the emissions from unsaturated wetlands.

Seasonal cycle
The mismatch between the region's carbon fluxes and saturated soil is not only spatial in nature; it is also temporal.The seasonal cycles of the meteorological forcings, wetland zone area fractions, and carbon fluxes, for the entire domain and the southern and northern halves are plotted in Fig. 8.While air temperature (T air , Fig. 7a-c) peaked in July, A sat (Fig. 7g-i) peaked in May-June, in response to snowmelt inputs and the drawdown from evapotranspiration (Melt and ET, respectively; Fig. 7d-f).Across the entire WSL and particularly in the south, all carbon fluxes (Fig. 7j-o) with the exception of

Conclusions References
Tables Figures

Back Close
Full CH 4 from the saturated wetlands (blue line, Fig. 7m-o) peaked in July (or July-August in the case of CH 4 ).In the south, where the bulk of the carbon fluxes were generated, July-August A sat values averaged only 0.2, or 30-40 % of their peak values.To assess the degree to which the July peaks in carbon fluxes were caused by inhibition of NPP and Rh during the June peak of A sat , we compared the seasonal cycles from the primary (distributed water table) simulations to those of control runs using a uniform water table scheme (so that A sat was always 0).The resulting carbon fluxes (denoted by dashed lines in Fig. 7j-o) were similar to the seasonal cycles of the noncontrol fluxes.Thus, it would appear that the carbon fluxes responded primarily to the July peak in T air , with A sat losing most of its potential influence over the fluxes by the time of their peak values.

Interannual variations (historic reconstruction)
As a consequence of the spatio-temporal mismatch between A sat and carbon fluxes, A sat had little influence over the interannual variability of carbon fluxes in the WSL.Correlations among the domain's annual carbon fluxes and JJA T air , precipitation, and A sat are listed in Table 5.As we might expect, JJA A sat displayed a negative correlation with JJA T air (via evaporative losses) of −0.63 and a strong positive correlation of 0.77 with JJA precipitation.Despite the fact that NPP and Rh both showed strong negative correlations with JJA A sat (−0.76 and −0.78, respectively), the net carbon flux (in which CH 4 plays only a minor role) showed very little correlation (0.16) with JJA A sat .Similarly, while CH 4 emissions from saturated and unsaturated wetlands displayed clear correlations with JJA A sat (0.72 and −0.49, respectively), their opposing signs and the larger area of unsaturated wetlands resulted in a very low correlation (−0.09) of total CH 4 with JJA A sat .GHWP also displayed little correlation (0.09) with JJA A sat , once again due to the opposing influences of its component fluxes.As a result, the domain's carbon fluxes displayed stronger correlations with JJA T air than with either JJA A sat or JJA precipitation.Control runs using a uniform water table scheme showed similar correlations between carbon fluxes and JJA T air and precipitation.

Conclusions References
Tables Figures

Back Close
Full Historical time series of these terms, plotted in Fig. 8, tell a similar story.JJA T air (Fig. 8a) displayed a general upward trend between the mid-1960s and 2010; performing a Mann-Kendall trend test on the time series for all segments of at least 20 yr in length yielded positive trends at the 95 % confidence level ranging from 0.023 C yr −1 to 0.080 C yr −1 for most segments having a starting year in the 1960s and an ending year between 1990 and 2010.Fewer trends in precipitation, evapotranspiration, and snowmelt (Fig. 8b) were significant over this period, but precipitation did have negative trends ranging from −1.1 to −2.2 mm yr −1 for some segments beginning in the 1960s and ending in the 1990s.JJA A sat (Fig. 8c) had significant negative trends ranging from −0.0004 to −0.0018 yr −1 for roughly half of the segments starting in the 1960s and ending between 1990 and 2010.NPP and Rh (Fig. 8d) both had significant positive trends ranging from 0.55 to 1.7 Tg C yr −1 and 0.47 to 1.2 Tg C yr −1 , respectively, over most segments in this period.Because control runs with a uniform water table also yielded similar positive trends in NPP and Rh (Fig. 8d, dashed lines), these trends were presumably in response to the positive trends in T air .Yet NEE (Fig. 8e) did not have any significant trends in segments starting in the 1960s, although NEE did have significant trends ranging from −0.09 to −0.23 Tg C yr −1 for some segments starting in the 1950s and ending in the 1990s.Similarly, CH 4 emissions from saturated and unsaturated wetlands (Fig. 8f) had significant trends of −0.005 to −0.023 Tg CH 4 yr −1 and 0.009 to 0.017 Tg CH 4 yr −1 , respectively, for many segments starting in the 1960s, presumably in response to the decline in JJA A sat , but total CH 4 emissions (Fig. 8f) and GHWP (Fig. 8g) do not had any significant trends.

Discussion
The most striking result of our work is the markedly different spatial distribution of methane emissions in the WSL given by this study, in comparison with most previous to North in the WSL, with the most pronounced gradient occurring between 58 • and 62 • N.This region coincides with both the southern limit of permafrost and a zone of large fractional extents of lakes, inundation, and saturated wetlands in the center of the WSL.As shown in Figs.7d-f, we were unable to match this drop in methane emissions without both (a) using a less productive methane emissions parameter set in the permafrost zone and (b) accounting for emissions from unsaturated wetlands.These factors can plausibly explain the differences between our spatial distribution and those of most previous estimates, which did not have access to the observations of Glagolev et al. (2011).For example, most previous bottom-up estimates (e.g.Schuldt et al., 2012;Zhu et al., 2012;Ito and Inatomi, 2012;Meng et al., 2012;Spahni et al., 2011;Riley et al., 2011;Ringeval et al., 2011Ringeval et al., , 2010;;Petrescu et al., 2010;Zhuang et al., 2004;Gedney et al., 2004;Fung et al., 1991) applied the same methane emissions parameters to all boreal wetlands (although Zhuang et al. and Zhu et al. distinguished between boreal and tundra wetlands), and consequently overestimated methane emissions in the central WSL.Curiously, Shindell et al. (2004), using the results of Walter et al. (2001), also overestimated methane emissions in the central WSL, despite allowing one of the methane emissions parameters (r0) to vary spatially as a function of mean annual temperature and NPP.The reason may be because this function was derived from a relatively sparse set of global in situ observations.In addition, some studies (e.g.Ringeval et al., 2010Ringeval et al., , 2011;;Gedney et al., 2004) employed wet-dry schemes that ignored the contributions from unsaturated wetlands (although Ringeval et al. did account for emissions from areas where the water table was within 5 cm of the surface), thereby underestimating emissions in the southern WSL.
Two previous estimates that did agree generally with our spatial distribution were those of Eliseev et al. (2008)  to a process-based model, which accounted for emissions from unsaturated wetlands, but applied a uniform methane emissions parameter set to all boreal wetlands.One of the five wetland area estimates did not overestimate methane emissions from the central WSL, but this seems to reflect the prescribed wetland areas rather than the model formulation.
This geographic discrepancy could have important consequences for predicting the response of wetland methane emissions in the WSL (and elsewhere in the high latitudes) to future climate change.Recently, concern has arisen in the cryosphere community that the thawing of permafrost could place previously frozen soil carbon at risk of decomposition (Koven et al., 2011;Schaefer et al., 2011;Walter et al., 2006).
How much of this carbon is respired anaerobically as methane depends on rates of methanogenesis and oxidation, estimates of which depend on assumed model parameters.Given that our single-parameter-set control simulation (Fig. 7d) yielded annual CH 4 emissions 2-3 times larger than our normal simulations in the central WSL (Fig. 7c) for the same amount of labile carbon, use of the same parameters in the permafrost zone as in the south could therefore lead to as much as a 2-to 3-fold overestimation of the response to permafrost thaw.While this assumes that methane emissions parameters (which depend to some extent on the composition of the soil microbial communities, which, in turn, may depend on environmental conditions) will remain constant in time, most studies to date have made such an assumption.
It is not clear how best to represent the spatial variation evident in the Glagolev et al. methane flux observations in a global process-based model.Our wetland LAI threshold of 2.5 for the more productive southern CH 4 parameter set was intended to contain the productive eutrophic mires found in the subtaiga, southern taiga, and middle taiga zones (Glagolev et al., 2011).At the very least, this suggests that global models may need to apply distinct methane emissions parameter sets at eutrophic mires, using specialized wetland land cover classifications similar to Peregon et al. (2009).A more process-based approach may require incorporating other spatially-varying factors, including soil pH, redox potential, substrate quality, and nutrient concentrations,

Conclusions References
Tables Figures

Back Close
Full although soil pH may play only a small role in spatial variability, due to the adaptations of microbial communities to local average pH (Glagolev, 2004).Some existing models already account for various combinations of some of these factors (e.g.Zhuang et al., 2004;Wania et al., 2010;Riley et al., 2011;Spahni et al., 2011;Meng et al., 2012;Zhu et al., 2012), but the inability of these models to reproduce the observed spatial distribution of emissions suggests that incorporating these factors may not be sufficient.Similarly, the low emissions under inundated conditions observed by  may indicate another process that is not currently accounted for in large-scale models and may help explain the spatial distribution of emissions.Carbon fluxes from our control simulations employing a uniform water table (for which A sat was always 0) differed only slightly from our fully-distributed simulations.However, it may be dangerous to interpret this as a validation of uniform water table schemes.The lack of influence of fractional saturation over current carbon fluxes appears to be the result of mismatches between the current spatial and temporal distributions of saturated soil and carbon fluxes.The current configurations of these spatial distributions indicate a substantial cumulative reaction to thousands of years of NPP inhibition under saturated conditions, which may not persist under projected future changes in the WSL' s climate.Indeed, Bohn and Lettenmaier (2010) found that the differences between the uniform and distributed water table schemes were most pronounced in their response to the climate projected for the end of the 21st century, rather than present-day fluxes.Nevertheless, it is important to constrain the fractional saturated area, if for no other reason than it is one of the few large-scale observations that can be used to calibrate a model's soil moisture storage and water table depth (to which carbon fluxes are extremely sensitive).However, as Fig. 3a shows, passive-active microwave products such as Schroeder et al. (2010) appear to be primarily sensitive to inundated extent, which is typically much smaller than the extent of saturated soil.dropped to 0 in many locations south of 60 • N. Thus, neither product alone is a reliable measure of saturated soil extent.Therefore, it would be very helpful to the modeling community for remote sensing specialists to generate products that combine both passive microwave and other (e.g.SAR-based) data.

Conclusions
We examined the role of spatial heterogeneity of surface and sub-surface water on the carbon budget of the wetlands of the WSL.We conclude that: -Sub-surface heterogeneity played an important role in both the overall magnitude and spatial distribution of estimates of the region's carbon fluxes, whereas surface heterogeneity had little overall effect.This was primarily because the bulk of the region's carbon fluxes occurred in the portion of the region where fractional saturated areas were lowest.
-To reproduce the observed spatial pattern of CH 4 emissions, in which very little CH 4 is emitted north of 60 • N, it was necessary to account for CH 4 emissions from unsaturated wetlands, and to use a methane model parameter set that reduced estimated CH 4 emissions in the northern half of the domain.
-Previous estimates of the response of the WSL to thawing permafrost may have overestimated future increases in methane emissions in the permafrost zone.

BGD Introduction
Full

A1 Hydrology
We used a modified version of VIC 4.1.2,an early version of which was described in Bohn et al. (2007).The model divides the land surface of each grid cell into separate land cover "tiles"; one of these tiles may be reserved to contain lakes and wetlands.Within the lake-wetland tile (with fractional area A wet ), lakes and wetlands are simulated as a continuous, connected system representing a single composite of all lakes and wetlands within the grid cell (Fig. 2).At any given time, some portion of the lakewetland system (A inund ) may be inundated with standing water on the surface, while the remainder of the system (A wet -A inund ) is exposed.The inundated portion may include both "permanent" lakes (A lake ) and seasonally-inundated wetlands (A inund -A lake ).
As described in Bowling and Lettenmaier (2010), A inund varies with time as a function of topography (or lake basin bathymetry) and the volume of impounded surface water.Surface runoff and sub-surface drainage from the exposed wetland flow into the inundated portion.Sub-surface drainage from the inundated portion of the lake-wetland flows into the channel network.All sub-surface drainage is controlled by VIC's baseflow parameters Ws, Ds, and Dsmax.When the water level in the inundated portion is above the lip of the permanent lake basin, the impounded water may also flow over the lip into the local channel network (controlled by the effective outlet width parameter, wfrac).If the inundated portion expands into the unsaturated wetland, some of the impounded water must recharge the newly-flooded wetland area until its soil is saturated.
Within the exposed wetland, field observations (e.g.Eppinga et al., 2008;Glagolev et al., 2011) indicate that microtopography exerts the single strongest control on local water table depth, with variations in water table depth over small scales nearly equal to variations in surface elevation relative to a datum.Therefore, we assume the landscape is composed of a mix of identical mounds (ridges) and depressions (hollows),

Conclusions References
Tables Figures

Back Close
Full as shown in Fig. 2. Following the observations of Eppinga et al. (2008), we assume the soil surface elevation within a hollow uniformly spans a range of 20 cm, while the surface of a ridge rises to a maximum to 50 cm above the edge of the hollow, for a total elevation range of 70 cm.The fraction of the landscape covered by ridges (f ridge ) is a calibration parameter.This topographic distribution is then sampled at regular intervals.Local water table elevation Z WT i is then computed at these same points via the equation where Z WT is the spatial mean water table elevation, Z surf i is the local microtopographic elevation, and Z surf is the spatial mean of the microtopographic elevation distribution.
The spatial mean water table depth is computed from total soil column water storage by following the method of Frolking et al. (2002), in which the water in peat soils is assumed to be in equilibrium between the forces of gravity and matric tension.
The locus of all points where the water table is at or above the surface is called the saturated zone, with fractional area A sat .Thus, A sat is the time-varying sum of permanent lake area, inundated wetlands, and saturated exposed wetlands: where f sat is the mean saturated area fraction across all ridge-hollow pairs in the exposed wetland.

A2 Biogeochemistry
In the modified VIC model, photosynthesis and aerobic soil respiration are simulated in the exposed wetland (and set to 0 under inundation Full where w min is 0, w max is 1, w opt is the optimal soil moisture content of 0.5 (at which peak respiration occurs), and r sat is the ratio of the respiration rate at saturation (w = 1) to the peak respiration rate (w = 0.5).k and r sat are calibration parameters.The total soil column Rh is integral of Eq. ( A5) across all 10 cm intervals between the surface and 2.5 m depth.The litter C pool exists only at the surface, while the intermediate and slow C pools co-exist throughout the remainder of the soil column.
Wetland methane emissions are computed using the model of Walter and Heimann (2000), driven by daily NPP, soil temperature profile, and water table depths from VIC (as described in Bohn et al., 2007 andBohn andLettenmaier, 2010).Methane emissions are computed separately for each point in the water table distribution.Lake methane emissions are not simulated.We modified the model of Walter and Heimann (2000) to respond to spatial variation in NPP.In the original model, the time series f in of substrate availability was given by: where NPP max is the local historical maximum daily NPP rate and f NPP (t) is a time series of substrate flux into the soil equal to daily NPP(t) during the growing season and a time-varying fraction of the previous growing season's total NPP during the subsequent winter.We replaced Eq. (A7) with:  Full  Full Methane flux observation sites (red circles) taken from Glagolev et al. (2011).Permafrost zones after Kremenetski et al. (2003).
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 5. As shown in Fig. 3c, for these optimal parameter values, VIC's A sat matched the PAL-SAR A sat reasonably well, with R 2 of 0.78 and bias −0.0015.The resulting JJA average simulated A sat values for the period 2001-2010, shown in Fig. 3d, tend to be much higher in the north of the WSL than in the south.Boxes in Fig. 3d delineate the grid cells used in the comparison of Fig. 3c.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 | 3 Results Discussion Paper | Discussion Paper | Discussion Paper | • N than to the north, due to the use of different CH 4 model parameter sets in the South and the North of the domain (addressed in Sect.3.3.).For C net (Fig. 5f), the largest carbon uptake rates (> 40 g C m −2 yr −1 ) occurred at the western edge of the domain, Discussion Paper | Discussion Paper | Discussion Paper | the two distributions compared favorably (R 2 = 0.53 for cells outside the region of the mismatch).The cause of the mismatch was not clear: the MODIS LAI values there were not substantially different from those of their neighbors, and none of our Monte Carlo simulations reproduced this feature.By design, our estimate of the total wetland carbon storage in the WSL, 70.7 Pg C (42.8-129 Pg C), matched the figure from Sheng et al. (2004) (70.2 Pg C) quite closely.3.3Spatial distribution of CH 4 emissions per grid cell areaFor comparison with our simulated CH 4 distribution, the spatial distributions ofGlagolev et al. (2011) andFung et al. (1991) are shown in Fig. 7a-b.The distribution of Glagolev et al., derived by mapping their observed CH 4 flux rates to the landcover classification of Peregon et al. (2009), places the vast majority of CH 4 emissions south of 60 • N, where wetlands are predominantly unsaturated.In contrast, the Fung et al. distribution (as well as many other subsequent estimates, e.g.Schuldt et al., 2013) places the majority of its emissions north of 60 • N, where saturated area fractions are larger and permafrost occurs.This difference in spatial distribution has important implications for the response of these wetlands to future climate change, as discussed in Sect. 4. To compare our simulated CH 4 emissions with these other estimates, we converted from flux per unit area of non-lake wetland to flux per unit area of the entire grid cell.The resulting distribution (Fig. 7c) placed the majority (70 %) of the domain's annual emissions south of 60 • N, with maximum values of 5-6 g CH 4 m −2 yr −1 confined to the Discussion Paper | Discussion Paper | Discussion Paper |region 65-85 • E, 57-59 • N. The broad features of this distribution agree with those estimated byGlagolev et al. (2011), although the emissions of Glagolev et al. were even more narrowly focused in the south than ours (likely due to our overestimation of saturated flux rates in the central WSL, shown in Fig.4h).Our underestimation of per-unit-area emissions in the far north (Fig.4f) did not produce large overall errors in the spatial distribution.Our estimate of total annual CH 4 emissions over the entire WSL, 3.65 Tg CH 4 yr −1 (1.69-5.96Tg CH 4 yr −1 ), is reasonably close to that ofGlagolev et al. (3.9  Tg CH 4 yr −1 ).This agreement is not surprising, since both studies used the same in situ observations from locations spanning the WSL.In contrast,Fung  et al. and Schuldt et al. estimated  substantially larger total annual emissions for the region (6-7 Tg CH 4 yr −1 ) To gain insight into the difference in spatial distribution between the Glagolev et al. and Fung et al. estimates, we performed control simulations in which either (a) the more productive CH 4 parameter set from the southern half of the domain was applied across the entire domain, or (b) the emissions from unsaturated wetlands were neglected (equivalent to a wet-dry scheme), or both (a) and (b) were applied.All Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | estimates.Both the in situ observations of CH 4 fluxes and the spatial distribution estimated by Glagolev et al. (2011) indicate a dramatic drop in emission rates from South BGD Discussion Paper | Discussion Paper | Discussion Paper | and one of the simulations of Petrescu et al. (2010).However, instead of accounting for emissions from unsaturated wetlands, Eliseev et al. assumed all wetlands were 100 % saturated, and prescribed wetland saturated extents statically from a map-based inventory rather than estimating them dynamically.Petrescu et al. (2010) used five different estimates of potential wetland areas as inputs Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | While the SSMI-based product of Papa et al. (2010) exhibited better agreement with the PALSAR saturated fraction in some portions of the central WSL, where forest cover is low, the Papa et al. product Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | w min )(w − w max ) (w − w min )(w − w max ) − (w − w opt ) 2 , w ≤ w opt w − w opt w max − w opt • r sat + (w − w min )(w − w max ) (w − w min )(w − w max ) − (w − w opt ) 2 , w > w opt (A6)

Table 1 .
Photosynthesis parameters for land cover classes in this study.Units are µmol CO 2 m −2 s −1 .

Table 3 .
Estimates of JJA average areas (km 2 ) of wetland zones, totaled across the WSL, over the period2001-2010.