Modeling micro-topographic controls on boreal peatland hydrology and methane fluxes

Small-scale surface heterogeneities can influence land-atmosphere fluxes and therefore carbon, water and energy budgets on a larger scale. This effect is of particular relevance for high-latitude ecosystems, because of the great amount of carbon stored in their soils. We introduce a novel micro-topographic model, the Hummock-Hollow (HH) model, which explicitly represents small-scale surface elevation changes. By computing the water table at the small scale, and by coupling the model with a process-based model for soil methane processes, we are able to model the effects of micro-topography on hydrology and methane emissions in a typical boreal peatland. In order to assess the effect of micro-topography on water the balance and methane emissions of the peatland we compare two versions of the model, one with a representation of micro-topography and a classical single-bucket model version, and show that the temporal variability in the model version with micro-topography performs better if compared with local data. Accounting for micro-topography almost triples the cumulative methane flux over the simulated time-slice. We found that the singlebucket model underestimates methane emissions because of its poor performance in representing hydrological dynamics. The HH model with micro-topography captures the spatial dynamics of water and methane fluxes, being able to identify the hotspots for methane emissions. The model also identifies a critical scale (0.01 km) which marks the minimal resolution for the explicit representation of micro-topography in larger-scale models.


Introduction
Peatlands cover only about 3 % of the global land surface (Wieder et al., 2009), but they play a fundamental role in the global carbon cycle (Blodau, 2002;Limpens et al., 2008).In boreal latitudes peatlands and wetlands are one of the major natural sources of methane to the atmosphere (e.g., Bousquet et al., 2006).During the Holocene peatlands have functioned as a sink of atmospheric carbon (Smith et al., 2004;Yu, 2012), Yu et al. (2011) recently estimated the amount of carbon stored in northern peatlands of about 547 (473-621) Pg, significantly larger than previous estimates of 270-370 Pg (e.g., Turunen et al., 2002).Recent efforts have tried to reproduce peatland and wetland extent and carbon accumulation in various Dynamic Global Vegetation Models (DGVMs), (i.e., Schuldt et al., 2013;Kleinen et al., 2012;Wania et al., 2009a, b).The WETland and Wetland CH 4 Inter-comparison of Models Project (WETCHIMP) (Bohn et al., 2015;Melton et al., 2013;Wania et al., 2013) revealed the large variability among the different models in wetland extent and in the parameterization of hydrological and biogeochemical processes such as methane emissions.It is also clear that all of these models lack the explicit representation of fine-scale heterogeneities and sub-grid processes.
We propose a novel method that takes into account subgrid scale processes, and directly assess their impact on peatland ecosystems, from greenhouse gas emissions, to the water budget.Previous studies have suggested that micro-F.Cresto Aleina et al.: Modeling micro-topographic controls topography in peatlands (micro-relief with a typical spatial scale of 1 m × 1 m) may influence GHG emissions (Baird et al., 2009b), but the extent of the micro-topography influence in land-atmosphere fluxes is yet to be determined.Small-scale processes can have significant effects in the peatland carbon cycle (Baird et al., 2009a), and local surface models (e.g., Nungesser, 2003) also highlighted how local surface heterogeneities matter for the water balance in northern peatlands.Acharya et al. (2015) recently linked the selforganization occurring in patterned peatlands to local smallscale interaction among micro-topography units.Other studies (Bohn and Lettenmaier, 2010;Bohn et al., 2007) examined the influence of water table heterogeneities on methane emissions using a TOPMODEL approach, which uses a modified water table in a single bucket grid cell (Beven and Kirkby, 1979).Bohn et al. (2013) and Wania et al. (2010) introduced a more physically based representation of hummocks and hollows, but none of these studies investigated the process-based representation of the lateral flow between hummocks and hollows.Cresto Aleina et al. (2013) instead showed that the importance of small-scale surface heterogeneities in estimating water table change in permafrostgenerated soil patterns.Observations in northern peatlands also showed the position of the water table has a fundamental control on greenhouse gas emissions, since it changes the depth of the oxic zone, i.e., the region where methane gas diffusing from below can be oxidized and therefore released as CO 2 instead (Couwenberg and Fritz, 2012).Because of the high global warming potential of CH 4 relative to CO 2 a robust estimation of methane emissions is essential to evaluate the climate impact of natural wetlands and peatlands (Kirschke et al., 2013).To compute a consistent greenhouse balance over the region, one should consider the small-scale properties and how the water table and the soil surface heterogeneously change within the environment (Bellisario et al., 1999;Camill and Clark, 1998;Law et al., 2002).Process-based models recently suggested that hydrological heterogeneities at the landscape scale between different wetland types (i.e., between fens and bogs) control different water table responses under a changing climate forcing (Gong et al., 2012).This phenomenon can potentially influence the carbon fluxes from peatlands at regional scales (Gong et al., 2013).Van der Ploeg et al. (2012) showed how micro-topography exerts dominant control in hydrological processes in wetlands.On the other hand, there has not yet been a quantification of the micro-topography effects on methane emissions, nor is there a proper way to represent these effects on larger-scale models.The general issue of scale interactions in the climate-biosphere system is therefore of particular interest in northern peatlands, where large emissions of greenhouse gases are influenced by the smallscale surface heterogeneities.
Here, we developed a mechanistic model operating at the landscape scale for a typical boreal peatland, in order to assess the importance of surface micro-topography for the wa-ter balance and the methane fluxes.These small-scale surface heterogeneities are typical in peatlands, and consist of elevated and relatively drier zones, called hummocks, and lower and relatively wetter zones, called hollows.We calibrated this landscape-scale model (Hummock-Hollow model, HH model hereafter) with data from an elevated bog in the Ust-Pojeg mire complex, in the Komi Republic, Russia.A number of recent studies have analyzed this site's peat characteristics and depth (Pluchon et al., 2014), and provide measures of fluxes of water vapor (Runkle et al., 2012), carbon dioxide (Schneider et al., 2012), andmethane (Gažovič et al., 2010;Wolf, 2009), as well as the energy and water balance (Runkle et al., 2014) and spatial distribution of dissolved organic carbon (DOC) (Avagyan et al., 2014b) and biogeochemical elements (Avagyan et al., 2014a).
In this paper we present the new model for peatland microtopography and the measurements we used for tuning and evaluating model performances.We then proceed to analyze the influence of micro-topography on water balance and methane processes comparing the output from the model including a representation of micro-topography with the output of a version of the model without any small-scale information.

Methods
We developed the Hummock-Hollow (HH) model to evaluate peatland micro-topographic controls on land-atmosphere fluxes.We represented the surface sub-grid scale heterogeneities with a square lattice.The model works at a landscape-scale of 1 km × 1 km, and each model grid cell represents a micro-topographic feature, namely a hummock or a hollow, with dimensions of 1 m × 1 m.We compute for each grid cell both water table balance and methane emissions.We explicitly represent the micro-topography characterizing the peatland surface, and therefore we parameterize the heterogeneous hydrological properties of the peatland soil (i.e., surface and the subsurface water fluxes).Due to such a fine representation of the micro-scale, we can upscale emissions and water balance at the landscape-scale by averaging the local quantities over the whole domain.To investigate the micro-topographic effect, we created a second version of the model (hereafter Single Bucket), in which we represent the whole peatland in a single bucket with parameters averaged over the whole 1 km × 1 km domain.

The HH Model
We tuned the model with micro-topographic data from a peatland of Northwest Russia, the Ust-Pojeg mire in the Komi Republic (61 • 56 N, 50 • 13 E, 119 m a.s.l.).Many studies have focused on this study site as a typical boreal peatland, as the mire complex displays different kinds of peatland types, from an ombrogenous bog, to a mineroge-nous fen, to a transitional zone with the surrounding forest.In order to simulate small-scale surface heterogeneities, we use field data to initialize the surface elevation in the model.We focus our modeling framework on the bog, in order to exclude the influence of the subsurface water input that characterizes the fen ecosystem.
We distinguish between two different elevations, a microtopographic one (i.e., the difference in elevation with respect to a surface level between hummocks and hollows), and a macro-topographic one, that takes into account the differences among grid cells due to the slope of the peatland.
We initialize the soil surface elevation with elevation data collected through surveying with a theodolite.The circumference and points on the height of each hummock and within representative hollows in a 40 m × 60 m grid were surveyed and a surface between them interpolated.In order to consistently represent the surface elevation, we inferred the statistical distribution of the field data.We fitted the histogram of elevational data (x) with different distributions (normal, logarithmic, 2-parameters gamma), and we found the best fit for a generalized 3-parameter gamma distribution multiplied by the maximum number of counts in the histogram: Values of parameters for Eq. ( 1) are reported in the following section.We randomly pick a value from the statistical distribution, and we assign it to each grid cell.We assume that with this procedure we statistically capture the peatland microtopography.If the grid cell at the position i, j has an elevation H i,j (height in Fig. 1) above a surface level H 0 = 20 cm then we assume that such a grid cell is a Hummock, otherwise it is a Hollow.
Along with the micro-topography, we initialize each grid cell with two other properties which affect methane production and hydrology respectively: peat depth pd(i, j ) and absolute elevation sl(i, j ).We model the peat profile (peat depth) according to in situ measurements, and we use this information in the methane emission model as a surrogate for the amount of carbon available for decomposition and GHG emissions.The equations used for pd initialization can be found in the Appendix.The second term sl(i, j ) parameterizes the macro-topographic elevation of each grid call with respect to the lowest part of the bog.Following observations, we assume the slope of the bog to have a uniform linear dependence on y, which represents the distance in meters from the origin.We therefore assign to each cell an absolute elevation according to the formula: where m y = 1000 m is the longitudinal dimension of the landscape-box representing our bog and HS = 3 m is the maximal height of the bog.The slope regulates water flow over the peatland (see Sect. 2.2.1).The total elevation of each grid cell is therefore the sum sl(i, j ) + H i,j .It is therefore able to resolve the micro-topographical features such as hummocks and hollows.Each grid cell has an elevation which is randomly assigned from the distribution of elevation data collected in the Ust-Pojeg mire complex in the Komi Republic, Russia.For each grid cell we simulate a dynamical water table, which changes with precipitation (P in the figure), evapotranspiration (ET), and lateral runoff (R).
In the Single Bucket model version all quantities are averaged, in order to assess the effect of the representation of micro-topography if compared to a mean field approximation.

Water table dynamics
We compute the water table position with respect to the surface level at each time step.The simulations start at the end of April and end at the end of October, when the cold season starts.After the snowmelt (not simulated) most of the landscape becomes inundated.For each cell at position i, j we compute the water balance with the following equation: where W i,j is the water table level in the grid cell at the position (i, j ) relative to the surface level, P is the precipitation input, ET is the evapotranspiration, R i,j is the lateral runoff, s i,j is the drainable porosity, and t is time, and the time step for the simulation is δt = 1 day.We randomly initialize water table level at the beginning of the simulations using values measured by Schneider et al. (2012) typical for hummocks and hollows: Drainable porosity in hummocks 0.8 - Drainable porosity in hollows 0.5 - Coefficient for Darcy's Law in hummocks 0.005 m d −1 Coefficient for Darcy's Law in hollows 0.01 m d −1 Coefficient for Manning's flow in hummocks 0.075 - Coefficient for Manning's flow in hollows 0.05surface level, s i,j = 1, whereas if the water table is below the surface, we assume it to be heterogeneous within the environment, and we choose two different values for hummocks and hollows.Values of s i,j are described in Table 1.By dynamically representing the water table at the micro-topographic level, we are able to distinguish the subsurface water flow in hummocks and hollows.We expect the spatial pattern of water table depth to change over the season, as the water flow, along with the evapotranspiration, makes the system progressively drier.

Lateral runoff
The lateral runoff term R is the computed water flux among the different grid cells, and it allows the exchange of water through the soil.We define it as R in i,j − R out i,j , where R in/out i,j is the sum of Darcy's and Manning's flows D and M: These terms represent the subsurface water flow and the overland water flow respectively and are computed in mm/day.We parameterize the subsurface water flow as The term m x/y i,j is the distance between the centers of cells i and j , and W i,j is the difference between the water table in the two adjacent grid cells.With respect to the classical expression for Darcy's Law, we then consider the term d i,j to represent the hydraulic conductivity, and its dimensions are length per time.It is spatially dependent, since we assume a different value for hummocks and hollows, as reported in Table 1.Parameters are chosen in the order of magnitude of the ones measured by Clymo (2004).Manning's formula describes instead the velocity of an overland flow driven by gravity.Following Manning (1891), we parameterize the surface flow as The overland water flow is also dependent on the soil heterogeneity because we assume a difference in the hummockhollow surface roughness due to the increased vascular plant cover in hummocks.In the equation, this term is represented by the dimensionless number Mc i,j , as displayed in Table 1.
In Eq. ( 7), sl 1 2 i,j represents the squared root of the difference in slope between the cells i and j while k is a conversion factor of 1 m 1/3 s −1 .The term Rad i,j is the hydraulic radius, defined as Rad i,j = A i,j /p i,j , where A i,j is the cross sectional area of flow, in our model: and p i,j is the wetted perimeter, or the perimeter of the cross sectional area A i,j at contact with water, in our model: where W i,j − H i,j is the elevation of water table with respect to the surface and m x/y i,j is the lateral extent of the grid cell.Parameters are chosen in the order of magnitude of the ones measured by Phillips and Tadayon (2006) for light and medium to dense shrubs, in the absence of specific values for hummocks and hollows.
Manning's flow occurs only if the water table is above the surface level, whereas Darcy's flow is continuous.This hydrological representation is one of the main differences between this new approach and the classical bucket model, and an important driver of the model's seasonal hydrological dynamics.The heterogenous surface and the interactions among the different grid cells represent at a fine spatial scale the interactions among hummocks and hollows in typical peatlands.Values for parameters are displayed in Table 1.
In the Single bucket configuration we used the same formulas, but with one single grid cell, i.e., i = j = 1.

Coupling with a process-based model for CH 4 emissions
We propose an explicit parameterization of methane fluxes, by coupling the micro-topographic water balance model to a process-based model for methane emissions, in order to more consistently quantify the effect of surface heterogeneities on GHG fluxes.The model developed by Walter and Heimann (2000) is quite a general model for methane emissions, and can be applied to peatlands in different environments.In particular, it is the same model which is built into some DGVMs (i.e., Kleinen et al., 2012).We tune the Walter and Heimann (2000) model parameters to perform in a typical peatland at the latitude of the Ust-Pojeg mire complex.The tuning parameter R 0 in the methane emission model has been adjusted depending on climatic condition.We set it to 0.30, following references in Walter and Heimann (2000).We couple the methane model at each grid cell, and we compute methane fluxes for each hummock and each hollow.We average over the whole landscape in order to upscale the local fluxes at the landscape-scale.
where W i,j is the water table depth with respect to the surface computed at each position i, j , z soil i,j is the soil depth, and NPP and T are at the daily timescale.The soil depth takes into account that each grid cell has a different peat depth.We sum the peat depth to the height of the acrotelm (the part of peat containing living plants, which we assume to be H i,j ).This quantity z soil i,j is a proxy for the amount of carbon available by adding to the peat depth the micro-topography height: In this application of the model we considered no ecological differences between hummocks and hollows, and therefore ignored potential differences in vegetation controls on methane emissions.

Sensitivity analysis
In order to test the robustness of our results, we change the key parameters of the model.We call hereafter HH Standard configuration the version of the model with parameter values described above.
In particular, we test our model by varying the grid size.In the HH Standard configuration the number of grid cells N = 10 6 , and the grid cells have a length size of m x/y = 1 m.We compare the output of the model with micro-topography with a classical Single Bucket model, i.e., a model where N = 1 and m x/y = 10 3 m.To test the scale dependency of the model, we perform simulations that gradually increase the number of grid cells (i.e., decreasing the grid cell size), and we compute the cumulative methane emissions over the whole season, from the end of April to the end of October.We aim to find a critical scale at which the model performance does not change qualitatively.
Another parameter to be tested is the drainable porosity s i,j , which plays a key role in the water balance, according to Eq. (3).By changing this value, the water table responds differently.In particular, for small values of the parameter, the water table variations within the soil are amplified.We also explore the importance of the difference in drainable porosity among hummock and hollows.
The peatland slope sl(i, j ) described in section 2.1 controls the surface flow, and in order to distinguish the effects of lateral flow and the ones of micro-topography we performed a further control simulation for sl = 0.
We also studied the dependence of the hydrological properties of the model by varying to the parameters in Table 1, which mainly influence the velocity of the surface and subsurface flow among the grid cells.We also ran the HH model forced by half the NPP values in the HH Standard configuration, and then we forced with NPP values which are twice the original ones in order to assess the goodness of our assumption on NPP values (see following sub-section).

Forcing data
The HH model is forced with prescribed precipitation and evapotranspiration (Eq.3).We use precipitation from Schneider et al. (2012), and evapotranspiration computed by Runkle et al. (2014).ET values in this paper were computed for the fen part of this mire based on eddy covariance measurements, but we use them directly for the bog part of the peatland.The variables are applied uniformly through the peatland, because we assume the spatial variation to be neglectable.
We force the model developed by Walter and Heimann (2000) with a prescribed time series of net primary productivity (NPP) and temperature (Eq.10).The time series for NPP are computed from simulations performed for the CMIP5 experiments with the MPI-ESM model at T63 resolution.We extracted NPP of C3 grasses for the grid cell which corresponds to the Ust-Pojeg mire.We are aware that moss has a different NPP than the typical plant functional type (PFT) representing C3 grasses, which is used in JSBACH, the land surface scheme of the MPI-ESM model.However, we use C3 grasses for simplicity, since mosses are not explicitly simulated in JSBACH.This approximation may introduce biases in estimations of methane emissions.We investigate the effect of this approximation in a sensitivity analysis (as described in the previous sub-section).In order to test this assumption, we investigated the sensitivity of out results by changing the NPP input values of Eq. ( 9).Potential changes in thermal insulation or carbon uptake due to moss coverage such as the ones highlighted by Porada et al. (2013) are beyond the purpose of this paper.
We extract daytime mean temperature from measurements by Runkle et al. (2014) in order to compare model output with observations.Further investigations of the potential influence of micro-topography on the energy balance and the land-atmosphere heat fluxes will address this potential source of differences between model and field measurements.

Model-data comparison
CH 4 fluxes were measured once a week from 27 April to 31 October 2008 applying a closed chamber approach (chamber dimensions: base 60 cm × 60 cm, height 25 cm).A total of 18 permanent measurement plots equipped with collars were established within the intensive study site in different microform types: two replicates each in ombrogenous hollows, lawns and hummocks, and three replicates each in minerogenous hollows, lawns and hummocks, and Carex rostrata lawns.The chamber was equipped with a fan to ensure an even mixing of the air inside the chamber and with a venting tube to avoid under-pressure during gas sampling and to allow the ambient pressure fluctuations to be trans- mitted into chamber headspace.Six air samples were taken from the chamber headspace during the 15-20 min chamber closure period using 60 mL plastic syringes.The air sample analysis was usually performed within the day following field-sampling with a gas chromatograph (GC, Hewlett Packard) equipped with a GFT PORAPAK a 80/100 (MESH-COND1900GC-015-9239, Hewlett Packard, USA) column and a flame-ionization detector (FID).Flux rates were calculated from the change of the CH 4 concentration in the chamber headspace over time by fitting a linear function by ordinary least-squares regression.

Results and discussion
The HH model allows us to study the potential landscapescale effects of micro-topography in a typical northern peatland.In the following sub-sections, we first discuss the statistics of the micro-topographic representation and we then present the effects of micro-topography on hydrology and methane emissions by comparing the novel approach to a classical bucket model and the performances of these two versions of the HH model against field measurements.

Micro-topography statistics
We compare the histogram of elevational data collected in the field with a 3-parameter gamma distribution (Eq.1).We confirm the goodness of the visual fit in Fig. 2 by running a Kolmogorov-Smirnoff test at a level of confidence of 95 %.The test shows that the population has no significant difference to the function f (x) in Eq. ( 1).We also tested other dis- The micro-topography affects water table position by delaying the runoff, because of the complex interactions among the grid cells, as resolved in a finer scale model.The shaded areas represent the standard deviation over 30 simulations (but their small size renders them invisible in the microtopography case).
tributions, such as normal, exponential, and lognormal, but none of them passed the Kolmogorov-Smirnoff test.We then assume that f (x) fits our data well enough to be used for our purposes as parameterization of the distribution for microtopographic elevation, and we proceed to initialize the microtopographic model by assigning at each grid cell a value randomly picked from f (x).The best fit parameters in Eq. ( 1) are (b = 5.8, q = 8.9, p = 1.5).It is worth noticing that for this peatland the distribution of surface elevation is not bimodal, as shown in previous studies (Eppinga et al., 2008).In fact, the two peaks in the histogram are too close to be resolved by the fitting distribution (Fig. 2).

Hydrology
Because of the random initialization we performed ensemble simulations with 30 ensemble members and we compared two versions of the model.Each ensemble member differs for initialized water table and, in the Microtopography version of the HH model, for micro-topography configuration.The ensemble members differ very little with each other in the Microtopography version, whereas the span of water table positions computed by the Single Bucket version is much larger.
The difference in water table dynamics between the two model versions is significant (Fig. 3).The shaded areas represent the standard deviations from the mean of 30 ensemble members and in the Microtopography version, this shaded area is so small that it is invisible in the Figure ( ±3 mm).This is due to the fact that the only random parameters are the micro-topographic height and the water table position at the beginning of the simulation.Since in the Microtopogra- With the HH model in the Microtopography version we are able to study the differences in water table position between hummocks and hollows, in order to compare the performance of the model against field data.In Fig. 4 we represent the distribution of simulated water table position relative to the surface level in box plots for hummocks (orange boxes), hollows (red boxes).To evaluate the effects of our hydrological model, we also plot the water table position computed by a realization of the HH model in the Single Bucket version (blue line).We compare model output against in situ measurements from the bog region (black crosses) for both micro-types.For completeness, we also compare the model performances against measurements from other regions of the Ust-Pojeg mire complex, i.e., the fen region (black circles in the figure).We chose these specific dates of Fig. 4 because of the relative abundance of chamber measurements of methane for both hummocks and hollows, which we will describe in the following sub-sections.The model performs generally well against the data, with most of the measurements among the whiskers of the box plots.The model seems to slightly overestimate the hummock water table towards the mid of August and to underestimate the hollow water table at the beginning of September, but the measurements still fall among the whiskers of the box.The model without microtopography, instead, consistently underestimates water table position.In particular, the computed water table is lower than the one of the deepest measurement for all analyzed dates.
In our simulation, the presence of outliers in the negative tail of the distribution shows that the data are skewed to the left.The skewness of both hummock and hollow distributions is negative for most of the simulation, becoming weakly positive for hummocks only at the beginning of the simulation, when water table is well above the surface level.Afterwards, the skewness of the hummock distribution reaches larger negative values than in the hollow distribution, as the hummocks are generally drier than the hollows.Accordingly, the excess kurtosis of the distributions in all cases is positive, since the tails of the histograms are fatter than they would be if data were normally distributed.In particular, the left tail becomes fatter as more and more grid cells become dry.

Methane emissions
By coupling the HH model with a process-based model for methane emissions, we are able to simulate the dynamics of methane over the warm season in the Ust-Pojeg mire complex.The HH model enables us also to distinguish among emissions from different micro-topographic units, in order to better compare the performances of the model against chamber measurements.
We compute methane emissions for each cell of the squared lattice.Due to the heterogeneous patterns of soil properties and water table position that we analyzed in the previous section, the emissions of methane are not uniform over the landscape.We study the impact of such heterogeneous emission pattern by distinguishing among hummock and hollow distributions of methane emissions.We compare the box plots of methane emissions at different time steps with chamber measurements taken in the bog part of the Ust-Pojeg mire complex (Wolf, 2009) put against in situ measurements from the bog region (black crosses) for both micro-types (Fig. 5).For completeness, we also compare the model performances against measurements from other regions of the Ust-Pojeg mire complex, i.e., the fen region (black circles in the figure).As expected, methane emissions differ between hummocks (orange boxes) and hollows (red boxes), being the latter ones generally wetter and thus displaying a shallower oxic zone.We also compare the averaged output of the model in the Single Bucket version (blue lines in the figure) against the distributions of methane emissions from the Microtopography version.In the data, the higher emissions from the wetter hollows in comparison to the drier hummocks are moderate evidence that hydrology is the main driver for the heterogeneity in methane emissions.
In this study we ignored potential chemical controls but the good agreement between data and model output suggests that our assumption is correct.The large increase in methane emissions at the end of July and at the beginning of August, as we can also see in the averaged fluxes (Fig. 6), is due to the combination of a higher temperature forcing and the presence of shallow oxic layer in most of the model grid cells in the Microtopography version.By the end of August, instead, despite the comparable water table position, temperatures are much lower.This phenomenon causes a significant decrease in methane emissions, only partially mirrored by chamber measurements.The HH model coupled with the Walter and Heimann (2000) methane emission model captures the general trend and the magnitude of methane emissions, but towards the end of the season it seems to fail in representing large methane emissions, in par-ticular from a hummock at the end of August and from two hollows in September, whereas in the rest of the season most of the measurements fall between, or near to, the whiskers of the boxes.Because of the generally good agreement of simulated water table depth with measurements, we exclude a bias coming from the hydrological model.The only other potential bias comes from the methane emission model, which seems to be overly sensitive to temperature variations.Despite such differences in fluxes, in general the measurements fall in the range of the simulated fluxes by the Microtopography version.The agreement is particularly good for hollows towards the center of the season, when methane emissions are higher.
The Single Bucket version represents an average over the whole region, but it nevertheless produces outputs which are not outside the range of the measurements.The inability of this version to represent a distribution of the large spread in methane emissions and in particular of hollows as methane emission hotspots, though, leads to an overall underestimation of the averaged fluxes (Fig. 6).The spatially explicit representation of the methane flux distribution, instead, is essential for the Microtopography version to better capture the magnitude of the measured fluxes in average.
The Microtopography version produces much larger fluxes than the Single Bucket version since, as we hypothesized, the latter version does not capture methane emission hotspots.The large peak toward the mid of July can be explained by the high temperatures coinciding with simulated wet conditions in the Microtopography version at the same time (Fig. 3).Because of the much drier average conditions in the Single Bucket version, the model is not able to capture such large spikes, which can be seen in the chamber measurements.Since the temperature forcing is the same as in the Microtopography version, the water table position is responsible for the poor performance of the Single Bucket version in representing the CH 4 fluxes.The water table is in this case much deeper than in the measurements, and this causes the methane produced to be partly oxidized and therefore the outgoing fluxes to be smaller than observed.The difference between the two model versions becomes even more striking by looking at the cumulative emissions over the whole warm season.The cumulative emissions in the two versions differ by almost a factor of 3, as the Single Bucket version produces (0.5424 ± 0.1931) × 10 4 mg m −2 of CH 4 , whereas the Microtopography version produces (1.5105 ± 0.069) × 10 4 mg m −2 of CH 4 .Micro-topography therefore controls methane emissions because of its influence on the peatland hydrology.

Critical scale and sensitivity to parameters
The results we presented were obtained with a specific choice of parameters.We tuned the model parameters with available values from the data set collected in the Ust-Pojeg mire complex and with standard values from the literature.In this sec- tion, we discuss the robustness of our results by changing the most important model parameters, starting with the grid cell size.The goal was to find a critical scale at which a finer representation of micro-topography did not significantly change the results.We therefore increased the number of grid cells from N = 1 (i.e., the model working in the Single Bucket version discussed in the previous sections with one grid cell of 1 km 2 size), to N = 10 6 (i.e., the model working in the Microtopography version, with a grid cell size of 1 m 2 ).In the previous section we chose the 1 m resolution in the Microtopography version because such resolution is approximatively the dimension of the micro-topographic features in the field.We computed then the cumulative CH 4 emissions for each simulation to test the dependence of emissions on grid cell size (Fig. 7).The cumulative emissions increase almost linearly for an increasing N , if N ≤ 10 2 .By increasing further the number of cells, the cumulative emissions stabilize after slightly decreasing for 10 2 < N ≤ 3 × 10 2 , but they do not largely change.Such different behavior depends on the different representation of water table dynamics, which delays the drying as the surface becomes more diverse in microtopography, i.e., if the number of grid cells increases.In particular, in the first 3 months of simulation the water is retained within the system as the number of hydrological interactions between the different grid cells increases, and the water table lays more and more in proximity of the average soil surface.This change explains the increase in emissions as N = 10 2 .By N > 10 2 , instead, the water table behavior asymptotically approaches water table dynamics for N = 10 6 , i.e., for model resolution of 1 m, hence changes in cumulative emissions are no longer significant.In our model, therefore, we can identify a critical scale for micro- N ≤ 10, i.e., decreasing cell size from 1 km 2 to 0.01 km 2 in a 1 km 2 domain.After this threshold, cumulative methane emissions stabilize as the number of grid cells increases.This phenomenon is mirrored by the average position of the water Among other parameters, the drainable porosity s i,j in Eq. ( 3) has a direct impact on water table position and therefore indirectly mon methane emissions.We chose the amount of water the peat soil can retain in the Standard configuration based on Kolka et al. (2011).This parameter can vary greatly in different peatlands, and therefore the values in our Standard configuration are an assumption we need to test.In order to assess the model sensitivity to this parameter we make drainable porosity uniform over the whole region, i.e., we do not distinguish between hummocks and hollows.We changed s i,j uniformly at steps of 0.1 from 0.2 to 1.0 to test how water table changes its position.The average position is generally lower with respect to the dynamics of Fig. 3, but the differences are not qualitatively significant, nor can we see any large effect in methane emissions.This is due to the methane emission model by Walter and Heimann (2000) being very sensitive to temperature variations, which in this test did not vary from the forcing of the Standard configuration.For very low values of s i,j , water table variation are larger and in general we observe both a deeper water table and lower methane emissions than in the Standard configuration, but the model still produces higher emissions if a representation of micro-topography is included.

Biogeosciences
The results proved to be robust also by changing the other parameters in Table 1 across a wide range of values.In particular, also for this analysis we eliminated the differences in hydrological parameters between hummocks and hollows, and the results were qualitatively the same as in the Standard configuration.
We also tested the influence of NPP forcing on model performance.In particular, we changed NPP values for 0.1 to 5 times the values in the Standard configuration.The differences that we saw in model output (see Appendix) were small or negligible and did not qualitatively alter the results.This finding confirms the robustness of our results, and the assumption that the bias introduced by not considering NPP for mosses, but only for C3 grasses was neglectable for the purposes of this study.
The control simulations for zero slope also confirmed the robustness of our results.The lateral flux and the water table decreases are less pronounced than in the Standard configuration case, but the differences between the HH model in the Single Bucket configuration and in the Microtopography are still significant (see Appendix).

Summary and Conclusions
We developed the Hummock-Hollow (HH) model, a new model representing the hydrology and the properties of the micro-topographic features typical of a boreal peatland, working at the resolution of 1 m 2 .This novel model presents a physical representation of the peatland micro-topography with the help of in situ measurements in the Ust-Pojeg mire complex in the Komi Republic, Russia.After inferring the statistical distribution of micro-topography data, we used this result to randomly assign elevation values at the grid cells.The explicit representation of the micro-topography allows the HH model to distinguish water table and methane fluxes among hummocks and hollows, thus identifying the role of the diverse micro-topographic features in water table and methane flux dynamics.To assess the effects of the microtopographic controls on these two observables, we created a model version which averages all quantities (Single Bucket), thus not distinguishing among hummocks and hollows, and which reproduces how a global or regional model would represent the landscape.We compared the output of the two model versions with in situ measurements of water table depth and methane fluxes.
Overall, the model version with micro-topographic representation performs better in comparison with hydrological data, as the water table position simulated by the Single Bucket version of the model is constantly deeper than the measurements.The Single Bucket version simulates a drier peatland, because the strong runoff washes away the water at the beginning of the simulation.The flow is instead diminished in the model version with a representation of microtopography, since the more rugged, hummocky surface delays the water discharge.This phenomenon allows for wetter conditions, leading to a general good agreement with field data in the Microtopography version working at 1 m 2 resolution.The HH model therefore correctly captures not only the averaged water table dynamic, but also the heterogeneity in water table depth distribution among hummocks and hollows, as micro-topography slows down the water runoff.
By changing the water table dynamics, the microtopography affects methane emissions.The water table position in respect to the surface level changes regulates the depth of the oxic zone, i.e., the region where methane can be oxidized and therefore the methane emissions are drastically reduced (as seen experimentally, e. g., in Couwenberg and Fritz, 2012).In our simulations, the Single Bucket version generally underestimates the averaged methane fluxes, because of the overly deep water table towards the central months of the simulations.The spatially explicit version of the HH model, instead, is able to produce an output of methane emission distributions, which as expected identifies the hollow grid cells as hotspots for methane emissions.
We progressively increased the scale of the model, i.e., reducing the number of grid cells, and we identified a critical scale at which the model results do not change for an increased resolution.This critical scale for the grid cells is 0.01 km 2 in the HH model, which is still too small for the investigated processes to be explicitly included in a global or regional model.Therefore we argue that further developments are needed towards a new parameterization of peatland surface which takes into account the upscaled effects of micro-topography.The identification of a critical scale for the representation of micro-topography on a global scale requires the application of the HH model to other peatlands, with a more structured and patterned micro-topography, which is the object of future investigations.
This last result limits the applicability of the model to landscape-scale studies, because of the computational nonfeasibility of including our findings directly into a global model.We tested the HH model only for one particular peatland, and even though we believe the peatland system we studied to represent a rather typical system in boreal peatlands, further research is needed to consistently assess the global relevance of our results.In particular, a necessary development to assess the larger scale relevance of our findings involves the application of the HH model to peatlands in different climatic regions.In this direction, the application of the HH model to a peatland with a more regularly patterned micro-topography, like the one described by Eppinga et al. (2008), can estimate the dependence of the micro-topography controls we described in the present paper on different microtopography configurations.
Further potential developments of the HH model include coupling the model to other process-based models for greenhouse gas emissions, CO 2 in particular, such as the one developed by Wania et al. (2010).Such a coupling would lead to the estimation of the micro-topography controls on the total carbon emissions of a typical boreal peatland.Moreover, our results are dependent on the particular choice of the process-based model for methane emission.The coupling with a model with an improved representation of litter chemistry could potentially have an influence not only on the quantitative results for methane emissions, but also on the definition of the critical scale.
Along this line of the development, the implementation of a peat accumulation module (e. g., the one presented by Kleinen et al., 2012) could potentially assess the microtopographic controls not only on hydrology and carbon emissions, but also on the long-term carbon cycle.Further model developments involve an explicit representation of the energy balance, in order to eliminate some of the biases introduced by the forcing time series, and to fully represent the dynamics of other processes we now ignore, i.e., microtopography controls on evapotranspiration, heat fluxes, and potential feedbacks among the different components of the model.The coupling with a model for moss dynamics (e.g., Porada et al., 2013) will make the model able to investigate dynamic and micro-topography dependent insulation properties of the soil.
The approach we developed for the HH model is not feasible to be directly applied in global modeling.The good agreement between the HH model and data shows that an explicit representation of micro-topography is fundamental in predicting landscape-scale hydrology and, as a secondary effect of surface heterogeneities, methane emissions.This study highlights the need of effective upscaling procedure in order to parameterize the effects of local surface heterogeneities across scales, ranging from the micro-scale (1 m) to the GCM scale (30-500 km).

Figure 1 .
Figure 1.Schematics of the HH model.The model represents a 1 × 1 km peatland, and works at a 1 m resolution.It is therefore able to resolve the micro-topographical features such as hummocks and hollows.Each grid cell has an elevation which is randomly assigned from the distribution of elevation data collected in the Ust-Pojeg mire complex in the Komi Republic, Russia.For each grid cell we simulate a dynamical water table, which changes with precipitation (P in the figure), evapotranspiration (ET), and lateral runoff (R).In the Single Bucket model version all quantities are averaged, in order to assess the effect of the representation of micro-topography if compared to a mean field approximation.

Figure 2 .
Figure 2. Comparison between the Ust-Pojeg mire topographic data collected in a field survey and a generalized three-parameters gamma distribution multiplied by maximum number of counts in the histogram.The good visual agreement is confirmed by the positive results from a Kolmogorov-Smirnoff test at 95% confidence (P ≥ 0.5).

Figure 3 .
Figure 3. Ensemble simulations of water table dynamics.The Microtopography version of the HH model (black line) is compared to the Single Bucket version (red line).The water table in the Microtopography version is averaged over the whole model region.The micro-topography affects water table position by delaying the runoff, because of the complex interactions among the grid cells, as resolved in a finer scale model.The shaded areas represent the standard deviation over 30 simulations (but their small size renders them invisible in the microtopography case).

Figure 6 .
Figure 6.Ensemble simulations of methane emissions averaged over the whole model domain.mWe compare the Microtopography version of the model (black line) with the Single Bucket model version (red line).We show the average of 30 ensemble runs (solid lines) and the shaded areas represent the standard deviation among the runs (but their small size renders them invisible in the microtopography case).

Figure 7 .
Figure 7. Dependence of the cumulative methane fluxes computed from the end of April to the end of October on the grid size of the micro-topographic model.The x axis shows the squared root of the number of cells N, as √ N = m x = m y , where m x and m y are the resolution horizontally and vertically, respectively.The x axis is on a logarithmic scale.Cumulative methane emissions greatly increase with increasing N if √N ≤ 10, i.e., decreasing cell size from 1 km 2 to 0.01 km 2 in a 1 km 2 domain.After this threshold, cumulative methane emissions stabilize as the number of grid cells increases.This phenomenon is mirrored by the average position of the water table dynamics (not shown), which by increasing N asymptotically approaches the water table dynamics for N = 10 6 .

Biogeosciences, 12, 5689-5704, 2015 www.biogeosciences.net/12/5689/2015/ F. Cresto Aleina et al.: Modeling micro-topographic controls 5695
Box plots of simulated water table depth in hummocks (orange boxes) and hollows (red boxes) against measurements from the bog (black crosses) and from the rest of the peatland (black circles).Grey crosses are outliers in the simulation.The blue line represents the output of the model without representation of microtopography.The box plot parameters regulating the length of the whiskers are the default values for covering 99.3 % of the data if the data are normally distributed.
phy version the model has a large number of cells (N = 10 6 ) all realizations are very similar to each other.Because of the initialization in both versions, the water table at the beginning of the simulations lies above the average surface level and floods the peatland.In the model version Single Bucket, the water flows out of the peatland much faster than in the Microtopography version.As a result, due to the strong surface (at the immediate beginning of the season) and subsurface flow, the water table drops quickly below 400 mm from the average surface level, thus increasing the oxic zone depth as the simulation proceeds.In the Microtopography version, instead, the sub-grid scale interactions delay the runoff.

net/12/5689/2015/ Biogeosciences, 12, 5689-5704, 2015 F. Cresto Aleina et al.: Modeling micro-topographic controls
Box plots of simulated methane emissions in hummocks (orange boxes) m and hollows (red boxes) against measurements from the bog (black crosses) and from the rest of the peatland (black circles).Grey crosses are outliers in the simulation.The blue line represents the output of the model without representation of microtopography.The box plot parameters regulating the length of the whiskers are the default values for covering 99.3 % of the data if the data are normally distributed . We compare model outwww.biogeosciences.