Modeling the large-scale effects of surface moisture heterogeneity on wetland carbon fluxes in the West Siberian Lowland

We used a process-based model to examine the role of spatial heterogeneity of surface and sub-surface water on the carbon budget of the wetlands of the West Siberian Lowland over the period 1948–2010. We found that, while surface heterogeneity (fractional saturated area) had little overall effect on estimates of the region’s carbon fluxes, subsurface heterogeneity (spatial variations in water table depth) played an important role in both the overall magnitude and spatial distribution of estimates of the region’s carbon fluxes. In particular, to reproduce the spatial pattern of CH 4 emissions recorded by intensive in situ observations across the domain, in which very little CH4 is emitted north of 60 ◦ N, it was necessary to (a) account for CH 4 emissions from unsaturated wetlands and (b) use spatially varying methane model parameters that reduced estimated CH 4 emissions in the northern (permafrost) half of the domain (and/or account for lower CH4 emissions under inundated conditions). Our results suggest that previous estimates of the response of these wetlands to thawing permafrost may have overestimated future increases in methane emissions in the permafrost zone.


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 historical 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 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.
The influence of soil moisture 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, the net impact of decomposition on climate.In northern peatlands, photosynthesis is also inhibited in some species (e.g., sphagnum) when the water table reaches the surface (Frolking et al., 2002).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) give 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 sub-surface 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 heterogeneous 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;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).

Modeling 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.2,as 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 (R h ), with inhibition of both of these fluxes under saturated conditions, and a heterogeneous water table scheme including 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 spin-up
Daily meteorological forcings for the period 1948-2010 were created by rescaling the NCAR/NCEP reanalysisderived fields of Sheffield et al. (2006) to match the monthly statistics of gridded monthly observations (precipitation and wind speed from Willmott and Matsuura, 2001;air Sheng et al. (2004).Methane flux observation sites (red circles) taken from Glagolev et al. (2011).Permafrost zones after Kremenetski et al. (2003).
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).
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 spin-up 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.One-sided leaf area index (LAI) values were prescribed by the monthly 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 R h 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/QuikSCATbased 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 Fig. 3b and d) from several dates in 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 of 0.0172, R 2 of 0.79).Therefore, we chose to treat the Schroeder et al. (2010) 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 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) digital elevation model (DEM) (Farr et al., 2007) south of 60 The remaining lake parameter w frac (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 June-July-August (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 (to which the water table depth is sensitive) needed to be determined: the fraction of the wetland covered by ridges (f ridge ) and the maximum sub-surface flow rate (D smax ).We held these parameters constant across the domain.We sampled f ridge and D smax 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 D smax , with the upper bound determined empirically by the maximum value of D smax 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 D smax equaled 0.5.As shown in Fig. 3c, for these optimal parameter values, VIC's A sat matched the PALSAR 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.

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 fac-tor 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 eddy covariance 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).Simulated soil temperature profiles (Fig. 4a-e) 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).Within each regional group, we divided the observations into four water table depth categories (Fig. 4f-j) and computed the mean of the log-transformed simulated and observed fluxes for each category.Our objective function was the joint likelihood of the categories' means across all groups under consideration.
We initially attempted to optimize a single set of methane emission parameter values across all groups (green symbols in Fig. 4f-j).This resulted in a pronounced negative bias in group 1 (Fig. 4e) and pronounced positive biases in groups 2 and 3 (Fig. 4f and g).Optimizing parameters separately for the northern (groups 1-3) and southern (groups 4-5) observations (blue symbols in Fig. 4f-j) remedied the positive biases in groups 2 and 3, with little change in the other groups.We suspect that our negative bias in group 1 might be due to a negative bias in our simulated NPP in the tundra.Peregon et al. (2008) found that tundra wetlands in the WSL (above 66 • N) exhibited annual NPP rates of approximately 1-5 times those of wetlands in the northern taiga (61-66 • N).In contrast, our simulated annual NPP rates in the tundra were, on average, 50 % of those in the northern taiga.Therefore, we multiplied our simulated NPP values at tundra sites by a factor of 3 and re-ran the single-and dual-parameter-set calibrations.The resulting single-and dual-parameter-set fluxes (yellow and gray symbols, respectively, in Fig. 4e-j) matched observed values in group 1 much better, with little change in other groups.We additionally tested alternate methods of computing mean annual soil temperature and allowing the efficiency of plant-aided transport to vary spatially, but effects on our results were small and are not shown here.Our simulated fluxes still failed to capture the substantial drop in methane fluxes between unsaturated and saturated conditions in groups 2-4; this will be addressed further in the discussion section.
Posterior distributions of the methane emission parameters for the various calibrations are listed in Table 3.A consistent difference between the single-and dual-parameter-set calibrations, with and without elevated tundra NPP values, lay in the values of rq10, which controls the sensitivity of methane production to seasonal changes in soil temperature (relative to annual average soil temperature).Calibrating the northern and southern sites separately, with or without elevated tundra NPP values, resulted in much lower rq10 values in the northern groups (2 to 5) than the southern groups (9 to 12).We will consider the implications of these differing parameter ranges in the dicussion section.
While the dual parameter set with elevated tundra NPP gave the best fit to observations, the cause of our low tundra NPP values was not clear (whether due to a low bias in unin-  [-] 0.00 0.25 0.50 k [-] 0.18 0.24 0.49 r sat [-] 0.15 0.34 0.50 hibited NPP, a high bias in inhibition under saturated conditions, or a high bias in inundated/saturated area).Therefore, we chose to use the dual parameter sets without modifying our simulated NPP values.We applied the northern parameter set to grid cells within the tundra and northern taiga, and the southern parameter set to all other grid cells.The boundary between these two sets of cells fell approximately at 61 • N (e.g., Glagolev et al., 2011) and coincided approximately with the southern limit of permafrost.

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 (via inhibition of NPP and Rh, and enhancement of CH 4 ), 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 4.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 km 2 and 332 000 km 2 , respectively).This represents 34 % (with 1st and 99th percentiles of 24 % and 47 %) of the non-lake wetland area of the WSL (and 13 % of the total area of the WSL).The remainder (66 %) (with 1st and 99th percentiles of 53 % and 76 %) of the non-lake wetland area was unsaturated.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 (non-lake) wetland carbon budget.As shown in Table 5, the total annual wetland NPP and R h 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 wetland methane flux of 3.6 Tg CH 4 yr −1 .Because NPP exceeds R h , 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 5); 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 total greenhouse warming potential (GHWP) of the region's wetlands 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 wetland NPP and R h 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 and wetland 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: b Area of wetland covered in standing water that is not permanent lake.c Area of wetland not covered in standing water but for which water table is at the soil surface.d Includes both "Non-lake Inundated" and "Exposed and Saturated".e Includes both "Non-lake Saturated" and "Permanent Lakes".
an increase of only 5 % in NPP, R h , and C net , and a decrease of only 2 % in CH 4 and GHWP, over the totals from the heterogeneous 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 wetland 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, R h , and soil carbon density (C dens ), shown in Fig. 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 with no seasonal inundation allowed (not plotted) yielded similar spatial distributions to the heterogeneous 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).However, 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 • N than to the north, due to the use of different CH 4 model pa- rameter sets in the south and the north of the domain (addressed in Sect.3.3).Even within the separate northern and southern halves of the domain, CH 4 exhibited only very weak positive correlations with A sat (0.27 and 0.18, respectively), due to the competing influence of NPP (which was negatively correlated with A sat ).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, 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 wetland GHWP (Fig. 5g) was more complex, with the large C net uptake rates in the west driving GHWP negative, and the large CH 4 emission 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 wetland soil carbon density (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 an undersimulation of carbon density in the tundra (66-70 • N, 75-90 • E), the two distributions compared favorably (R 2 = 0.53 for cells outside the region of the mismatch).The cause of the mismatch may have been an underestimation of tundra wetland NPP rates (see Sect.Fig. 6.Observed (from Sheng et al., 2004) and simulated soil carbon density, per unit wetland area.

Spatial distribution of CH 4 emissions per grid cell area
For comparison with our simulated CH 4 distribution, the spatial distributions of Glagolev et al. (2011) and Fung et al. (1991) are shown in Fig. 7a and b.The distribution of Glagolev et al. (2011) derived by mapping their observed CH 4 flux rates to the land cover classification of Peregon et al. (2008 and2009), places the vast majority of CH 4 emissions south of 60 • N, where wetlands are predominantly unsaturated.In contrast, the Fung et al. (1991) 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.Because these fluxes are from wetlands only, this distribution does not account for (typically small) methane oxidation in upland soils; therefore our estimated fluxes should be considered an upper bound on fluxes to the atmosphere.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 region 65-85 • E, 57-59 • N. The broad features of this distribution agree with those estimated by Glagolev et al. (2011), although the emissions of Glagolev et al. (2011) 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 tundra (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 of Glagolev et al. (2011) (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. (1991) and Schuldt et al. (2013) 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. (2011) and Fung et al. (1991) estimates, we performed control simulations in which either (a) the single CH 4 parameter set described in Sect.2.4 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 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 wetland 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. 8a-c) peaked in July, A sat (Fig. 8g-i) peaked in May-June, in response to snowmelt inputs and the drawdown from evapotranspiration (Melt and ET, respectively; Fig. 8d-f (Fig. 8j-o) with the exception of CH 4 from the saturated wetlands (blue line, Fig. 8m-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 R h during the June peak of A sat , we compared the seasonal cycles from the primary (heterogeneous 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. 8j-o) were similar to the seasonal cycles of the non-control 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 (historical 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 wetland 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 6.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 R h 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  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.Historical time series of these terms, plotted in Fig. 9, tell a similar story.JJA T air (Fig. 9a) 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. 9b) 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. 9c) 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 R h (Fig. 9d) 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 con-trol runs with a uniform water table also yielded similar positive trends in NPP and R h (Fig. 9d, dashed lines), these trends were presumably in response to the positive trends in T air .Yet NEE (Fig. 9e) 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. 9f) 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. 9f) and GHWP (Fig. 9g) did not have any significant trends.

Spatial variation of methane emission parameters
The most striking result of our work is the markedly different spatial distribution of methane emissions in the WSL given www.biogeosciences.net/10/6559/2013/Biogeosciences, 10, 6559-6576, 2013 by this study in comparison with most previous estimates.The spatial distribution estimated by Glagolev et al. (2011) indicates a dramatic drop in emission rates from south 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 widespread lakes, inundation, and saturated wetlands in the center of the WSL.This region and other areas along the southern limit of permafrost are of particular concern to the cryosphere community, because the thawing of permafrost (e.g., Smith et al., 2005) 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 will be respired anaerobically as methane depends on rates of methanogenesis and oxidation, estimates of which depend on assumed model processes and parameters.For example, using a single methane parameter set across the entire WSL (Fig. 7d) yielded annual methane emissions 2 to 3 times higher in the northern WSL than when using spatially varying methane parameters (Fig. 7c), for the same amount of input carbon substrate.
As shown in Fig. 7c-f, we were unable to reproduce the steep decline in methane emissions from south to north without both (a) using a less productive methane emission parameter set in the northern WSL (lowering northern emissions) and (b) accounting for emissions from unsaturated wetlands (raising southern emissions).These two 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., 2013;Zhu et al., 2012;Ito and Inatomi, 2012;Meng et al., 2012;Spahni et al., 2011;Riley et al., 2011;Ringeval et al., 2010Ringeval et al., , 2011;;Petrescu et al., 2010;Zhuang et al., 2004;Gedney et al., 2004;Fung et al., 1991) applied the same methane emission parameters to all boreal wetlands (although Zhuang et al. (2004) and Zhu et al. (2012) distinguished between boreal and tundra wetlands) and consequently overestimated methane emissions in the central WSL relative to our distribution and that of Glagolev et al. (2011).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 emission parameters (r0) to vary spatially as a function of mean annual temperature and NPP.The reason may be that 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 (Eliseev et al., 2008; and one of the simulations of Petrescu et al., 2010) did agree generally with our spatial distribution, but these appear to be the result of prescribed wetland extents rather than accounting for spatial variation in methane production or emissions from unsaturated wetlands.
Several possible underlying mechanisms could be responsible for the lower methane emissions in the northern WSL.One possibility is that the substantial differences we found between optimal rq10 values in the north and south represent real differences between the responses of soil microbial communities in those regions.This is somewhat corroborated by a recent synthesis of laboratory and field studies (Lupascu et al., 2012) that found higher mean Q 10 values for sphagnum-dominated wetlands than for sedge-dominated wetlands (means of 8 and 4.3, respectively), and higher (although not statistically significant) mean Q 10 values for nonpermafrost than permafrost wetlands.
Another possibility is that our lower rq10 values in the north were compensating for our model's failure to capture the relatively low methane fluxes observed under inundated conditions there.Similarly low methane emissions have been observed under inundated conditions in other studies.Strack et al. (2004) hypothesized that these low emissions were caused by oxidation in the water column above the soil surface, which our model did not account for.Oxidation in the water column therefore could play an important role in current and future methane emissions around the pan-Arctic, not only due to widespread inundation in permafrost lowlands such as the WSL and the Hudson Bay Lowland, but also under thermokarst conditions caused by thawing of permafrost (e.g., Smith et al., 2005).
A third possibility is that differences in emissions between the north and south are due to spatial variation in plant species assemblages.Several field studies (Glagolev et al., 2011;Levy et al., 2012;Olefeldt et al., 2013) have found strong links between plant species assemblages and wetland methane emissions.Similarly, modeling studies (e.g., Riley et al., 2011) have found that simulated methane emissions are very sensitive to species-specific vegetative parameters, such as the presence of aerenchyma, that affect both plantaided transport and root zone oxidation.
Other spatially varying factors, including soil pH, redox potential, substrate quality, and nutrient concentrations, could also play a role, 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.
Until the underlying mechanism is better understood, it is difficult to determine the applicability of our results outside the WSL.Despite having 750 observations, it is possible that sample error played a role in our parameter calibration, due to the highly heterogeneous nature (in time and space) of methane emissions.This is particularly likely in the tundra (represented by the approximately 25 observations, taken from a single wetland complex, comprising group 1).Comparisons with observations elsewhere in the high latitudes (e.g., Olefeldt et al., 2013;Levy et al., 2012) would be helpful in this regard.Similarly, 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.One approach could be to apply distinct methane emission parameter sets as a function of wetland type (e.g., eutrophic mires or ombrotrophic bogs) or plant species assemblage (e.g., sphagnum-dominated or sedge-dominated), using specialized wetland land cover classifications similar to Peregon et al. (2009).

Water table heterogeneity
Carbon fluxes from our control simulations employing a uniform water table (for which A sat was always 0) differed only slightly from our heterogeneous 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 both the dominance of unsaturated wetlands in the region's emissions and 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 heterogeneous 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 that 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.While the SSMI-based product of Papa et al. (2010) exhibited better agreement with the PAL-SAR saturated fraction in some portions of the central WSL, where forest cover is low, the Papa et al. product 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 the following: -Sub-surface moisture 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 spatially varying methane model parameters that reduced estimated CH 4 emissions in the northern half of the domain (and/or account for lower CH 4 emissions under inundated conditions).
-Previous estimates of the response of the WSL to thawing permafrost may have overestimated future increases in methane emissions in the permafrost zone.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 R h 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 f in (t) = 1 + f NPP (t)/NPP max , (A7) 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 timevarying fraction of the previous growing season's total NPP during the subsequent winter.We replaced Eq. (A7) with f in (t) = NPP max + NPP(t).(A8) As a consequence, the units of the tuning parameter r0 became (g C m −2 d −1 ) −1 .

Fig. 2 .
Fig. 2. Components of lake-wetland system considered in this study.

Fig. 3 .
Fig. 3. Observed and simulated fractional extents of inundation and saturation over the West Siberian Lowland.Units are the fraction of the entire grid cell area (wetland plus upland).

Fig. 4 .
Fig. 4. Observed and simulated soil temperature profiles (left column) and methane fluxes (right column) for each of the five observation groups delineated in Fig. 1 (rows 1-5, respectively)."Single" and "Dual" refer, respectively, to the use of a single methane model parameter set across the whole domain or separate parameter sets in the north and south."Tundra NPP" refers to simulations in which the simulated NPP of tundra sites was multiplied by a factor of 3.
2.4).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.

Fig. 7 .
Fig. 7. Spatial distributions of various estimates of annual methane emissions per unit area of grid cell (wetland plus upland).

Table 2 .
Posterior distributions of calibrated parameters for hydrology and soil respiration.

Table 3 .
Posterior distributions of calibrated methane emission parameters.
Walter and Heimann (2000)rom those inWalter and Heimann (2000); see Appendix A for details.

Table 4 .
Estimates of JJA average areas (km 2 ) of wetland zones, totaled across the WSL, over the period 2001-2010.
a Note: the areas of different categories within a given percentile do not, in general, sum to the total wetland area, because different categories cannot, in general, have the same percentile simultaneously.

Table 5 .
Annual average carbon budget terms over the period2001-2010.

Table 6 .
Temporal correlations of domain-wide annual carbon fluxes with JJA hydrologic conditions over the period 1948-2010.
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