A stochastic model for the polygonal tundra based on Poisson–Voronoi diagrams

. Subgrid processes occur in various ecosystems and landscapes but, because of their small scale, they are not represented or poorly parameterized in climate models. These local heterogeneities are often important or even fundamental for energy and carbon balances. This is especially true for northern peatlands and in particular for the polygonal tundra, where methane emissions are strongly inﬂuenced by spatial soil heterogeneities. We present a stochastic model for the surface topography of polygonal tundra using Poisson– Voronoi diagrams and we compare the results with available recent ﬁeld studies. We analyze seasonal dynamics of water table variations and the landscape response under different scenarios of precipitation income. We upscale methane ﬂuxes by using a simple idealized model for methane emission. Hydraulic interconnectivities and large-scale drainage may also be investigated through percolation properties and thresholds in the Voronoi graph. The model captures the main statistical characteristics of the landscape topography, such as polygon area and surface properties as well as the water balance. This approach enables us to statistically relate large-scale properties of the system to the main small-scale processes within the single polygons.


Introduction
Large-scale climate and land surface interactions are analyzed and predicted by global and regional models.One of the major issues these models have to face is the lack of representation of surface heterogeneities at subgrid scale (typically less than 100-50 km).How phenomena and mechanisms interact at different spatial scales is a challenging problem in climate science (see, e.g., the review of Rietkerk et al., 2011).In particular, the lack of cross-scale links in most general circulation models (GCMs) may lead to miss nonlinear feedbacks in climate-biogeosphere interactions, as highlighted by recent studies (Baudena et al., 2013;Janssen et al., 2008;Dekker et al., 2007;Scheffer et al., 2005;Pielke et al., 1998).This issue may cause a strong bias in climate models trying to compute accurate energy and carbon balances.This is particularly true for methane emissions in northern peatlands that contribute considerably to the Arctic carbon budget (Baird et al., 2009;Walter et al., 2006).In addition, thawing permafrost could lead to increased carbon decomposition and enhance methane emissions from these landscapes causing a positive feedback to climate change (O'Connor et al., 2010;Schuur et al., 2008;Christensen and Cox, 1995).
One important driver of methane emissions, for example, is water level.In polygonal tundra, the pronounced microtopography leads to a tessellate pattern of very contrasting water levels and hence strong small-scale variability of methane emissions.Previously, typical process-based wetland methane emission models (e.g., Walter et al., 1996;Walter and Heimann, 2000) were constructed based on plotscale process understanding and emission measurements but have subsequently been used for upscaling emissions to large areas using a mean water level.With water level often acting as a on-off switch for methane emissions (Christensen et al., 2001), the problem becomes obvious: where methane is emitted at the inundated microsites and not emitted at the dominating relatively drier sites, the mean water level in the model can either be near the surface (as in the inundated sites), leading to large emissions throughout the area, or it can be below the surface (as in the moist sites) leading to emissions close to zero.Both results would not represent the differentiated emission pattern of the polygonal tundra, also because water table position in respect to the ground surface is nonlinearly related to methane emissions.Zhang et al. (2012) also used a process-based model to upscale methane emissions from the plot to the landscape scale but ran the model separately for each land surface class as derived from high-resolution aerial imagery.This approach leads to a much better agreement with landscape-scale emissions as measured by eddy covariance than previous upscaling attempts with process-based models that ignored the spatial heterogeneity.The model representation of such heterogeneity is currently usually driven by the pixel size of available remote sensing products or the grid cell size of the respective model, but the real world functional or structural differences between individual components of a heterogeneous domain may actually be on a very different scale.Preserving their information content becomes particularly important where nonlinear relationships are involved in the represented processes -such as those driving greenhouse gas emissions (Stoy et al., 2009;Mohammed et al., 2012).
Arctic lowland landscapes underlain by permafrost are typically characterized by patterned ground such as icewedge polygonal tundra.They cover approximately 5-10 % of Earth's land surface, where they play a dominant role in determining surface morphology, drainage, and patterns of vegetation (French, 2007).Low-center polygonal tundra typically consists of elevated comparatively dry rims and lower centers that can become wet if the water table rises to the ground surface level or above.Thermally induced cracking in the soil during winters with large temperature drops generates the polygonal patterned ground.When the temperature rises and snow cover melts, water infiltrates into those cracks, and refreezes.Over the years, as this process is cyclically iterated, the terrain is deformed due to the formation of vertical masses of ice, called ice wedges.Aside from the wedges, terrain is elevated because it is pushed upwards, forming the elevated rims typical of this landscape.Incoming water (precipitation or snow melt) is then trapped by the rims in the polygon center, resulting in a water table level that dynamically responds to climate and weather.This environment is typical for high-latitude permafrost areas of Alaska, Canada, and Siberia.Ice-wedge polygon features have been observed also on Mars (Haltigin et al., 2011;Levy et al., 2010).The surface morphology of dry rims and wet centers determines drainage and vegetation patterns in the polygonal landscapes.They have therefore been the focus of extensive field studies on experimental sites on Samoylov Island in the Lena river delta in Siberia (an aerial picture of the study site is shown in Fig. 1).Previous works involved closed-chamber and eddy covariance measurements of methane fluxes (see e.g., Sachs et al., 2008Sachs et al., , 2010;;Wille et al., 2008;Kutzbach et al., 2004), as well as measurements and modeling of hydrological properties (Boike et al., 2008), and of surface energy balance (Langer et al., 2011a,b) in this typical periglacial landscape.Recent studies with the help of remote sensing data were able to capture and analyze surface heterogeneity (land cover and surface temperature) and its importance for the landatmosphere water fluxes (Muster et al., 2012;Langer et al., 2010).Theoretical efforts have also managed to quantitatively classify ice-wedge polygons and other permafrost patterns in Alaska using Minkowski density functions, quantifying binary patterns from aerial photographs (Roth et al., 2005).
In general, stochastic models may be useful tools to link different scales.Such an approach has been widely used in modeling physical, biological, and ecological phenomena, seeking to represent only the main processes and observable properties, and replacing complex dynamics with random processes (see e.g., Dieckmann et al., 2000, and references therein).To model large-scale (e.g., 100 km 2 ) seasonal dynamics of the small-scale-driven greenhouse gas emissions from such a landscape, it is impossible to limit the description to an accurate model for a single polygon (typical dimensions in the order of 10 m in diameter) due to the wide internal variability of the landscape.In particular, this mechanistic approach would require a large number of parameters and computational power in order to represent polygon variability in size, water table position with respect to the surface, and response to climatic forcing.In this work, we provide a stochastic method that takes into account such variability, using properties of Voronoi polygons.
Voronoi diagrams (or tessellations) have been applied to various fields, from astronomy to biology, forestry, and crystallography.For comprehensive reviews, see, e.g., Lucarini (2009), Okabe et al. (2000), and reference therein.Their wide range of applications derives from the fairly simple concept of partitioning the space behind them and from their adaptability.In ecology and Earth science, they have been used to describe ecological subdivision of territories (Hasegawa and Tanemura, 1976), hydrology and precipitation, geology, and recently also as a framework in order to integrate climatic variables over an irregular geographic region (Lucarini et al., 2008).

The model
We develop a stochastic model that is able to reproduce the main statistical characteristics of the low-center polygonal tundra (i.e., surface topography, area of different soil surfaces, water table height).With the present model we do not seek to explain surface pattern formation: our aim is, instead, to find a robust way to upscale land-atmosphere greenhouse gas fluxes and the landscape response to climatic forcings that are coupled to these statistics.We limit our description to the landscape's present-day state, and we represent neither polygon formation nor degradation.
We tessellate the plane with random polygons that are able to reproduce the main general features of ice-wedge polygons.The model represents, analyzes, and considers the statistical characteristics of a region covered by ice-wedge polygons and describes the landscape-scale (about 1 km 2 or more) response to climatic forcing.In order to represent icewedge polygonal tundra we use a random tessellation of the plane, so-called Poisson-Voronoi diagrams (Okabe et al., 2000).

Poisson-Voronoi diagrams
Let us populate a Euclidean plane with an at most countable number of distinct points.We then associate all regions in that space with the closest member(s) of the point set with respect to the Euclidean distance.The result is a set of regions that cover the whole plane, and each one of them is associated with a member of the generating points.This tessellation is called planar ordinary Voronoi diagram, generated by the point set.The regions constituting the Voronoi diagram are the Voronoi polygons.More rigorously, let P = {x 1 , . .., x n } be the generating point set, where 2 ≤ n ≤ ∞ and x i = x j for i = j , with i and j integers.The (ordinary) Voronoi polygon associated with x i is defined as and the set given by is the Voronoi tessellation.In our case, the generating process P is a homogeneous Poisson point process P (λ), where λ is the intensity of the process.Therefore, we call the resulting tessellation a Poisson-Voronoi diagram (PVD).Parameter λ basically controls the density of the Poisson process.On a square L × L, the intensity of the process defined as where E[n] is the mean of the number n of generated points, which follows a Poisson distribution.Parameter λ directly regulates the polygon sizes.We tuned this quantity with data from image analysis of the landscape (Muster et al., 2012).From Eq. (3), L = 1 km and λ = 7500.Since λL 2 1 we are in the thermodynamic limit and the boundary effects can be neglected.A model realization of PVD is shown in Fig. 2. The area of polygons in a Poisson-Voronoi diagram has been widely analyzed by numerous studies (Okabe et al., 2000).Direct numerical investigations showed that 3-parameter generalized gamma distributions fit quite well the statistical properties of the PVD (Miles and Maillardet, 1982): Other distributions have also been used in the literature to describe the PVD area distribution such as lognormal and 2parameter gamma (from Eq. 4, one obtains the 2-parameter gamma distribution if r = 1), as reported by Okabe et al. (2000) and reference therein.

Idealized polygons
Ice-wedge polygons are composed of two main regions: elevated drier rims and lower wet centers.In our model where A pol is the total polygon area, and r and C represent the area covered by elevated rims and the one covered by low centers, respectively.Parameter q is a random quantity varying between 0.35 and 0.7, which we tune according to observation data (Muster et al., 2012).Each Voronoi region represents an ice-wedge polygon and its vertical structure consists of a column, which parameterizes the water table level W t , the thaw depth TD, and the low polygon center surface level S. Water table is a dynamical variable and responds to climatic forcing, whereas thaw depth follows a prescribed behavior during the seasonal cycle (see Sect. 2.3 for the parameterization dynamics).Surface level, polygon depth, thaw depth, and water table level are initialized randomly: we tune their mean values with data from available field observations.All quantities are computed from the top elevation of rims, as shown in Fig. 3.
If the water table rises above the center surface, area C increases as well, following the equation where C is the polygon center area after the water table rise δW and α is the angle between elevated rims and polygon center surface, parameterized as which takes into account that the larger the rims, the steeper polygon walls become, as suggested by observations.Factor π 4 is needed not to exceed the limit α lim = π 2 .Computing water table (W t ) variations with respect to the polygon surface (S) is essential to estimate methane emissions (Sachs et al., 2008;Wille et al., 2008;Kutzbach et al., 2004).We distinguish between three kinds of surface types, depending on the position of the water table with respect to the polygon center surface.If where threshold = 10 cm has been inferred from observations of small-scale methane emissions (Sachs et al., 2010).They are controlled by the site's water table height, and different surface types correspond to different emitting characteristics.In particular, methane emissions from saturated centers are an order of magnitude larger than in the case of moist centers or than emissions from dry rims.

Water table dynamics
We simulate water focus on the effect of surface heterogeneity on GHG emissions, precipitation and evapotranspiration are parameterized as uniform over the whole area, i.e., we do not apply any downscaling of climatic forcing.In particular, the two variables follow a simple random parameterization.In this way, we are able to simulate both dry and wet summers.We also take into account that the water income in polygon centers does not only depend on precipitation income directly on the area C (Eq. 5) but some of the precipitation falling on the rims r also percolates through the soil and has an impact on the water balance.In particular, we assume that half of the precipitation falling on the rims percolates or flows to the polygon centers.We therefore consider an additional water income term.We generate a random number of rainy days (between 5 and 35).If the day t is a rainy day (r.d.), then where precipitation P is defined as where R p is a randomly generated number, which tunes the amplitude of precipitation events in the part of the summer season when main precipitation events are supposed to occur, and it is computed in mm day −1 .In particular, changing R p changes drastically the amount of incoming water.T = 30 days is hereinafter a time constant.Term S t represents the amount of water stored in the unsaturated terrain in the rims that does not contribute to water table variations at each time step.After large precipitation events, not all the precipitation flows from the rims to the water table centers but part of it fills the unsaturated terrain of the elevated rims.This water balance term has been observed in field campaigns and it is supposed to have a stronger influence at the beginning of our simulation time slice (July), when the upper layer of the soil is unsaturated, whereas storage capacity decreases with time until the whole terrain is saturated.In our model we parameterize storage as an exponential function: where t is time and τ = 45 days is the characteristic time at which S is reduced to (P −ET)/e.We assume precipitation and evapotranspiration to be uniform over the center and the rims.
We parameterize ET as where ET p follows We chose parameterization of ET and P in order to reproduce data from observation campaigns from Boike et al. (2008) and Langer et al. (2011b), which describe how evapotranspiration is linked to the landscape's the energy balance.
Thaw depth, on the other hand, follows a seasonal trend: its level depends on the time of the year the model is run at.With progress of the summer, the depth of the seasonally thawed layer (frost table depth) increases, reaching its maximum at the beginning of September.We keep thaw depth dynamics prescribed, i.e., we do not consider how temperature variations may affect the environment.From available observations, we assume thaw depth to behave as follows: During spring/early summer the ground is still frozen and water is kept inside the polygons.With progression of thaw depth during summer, neighboring polygons can become hydraulically connected.Field experiments show that water that was retained inside a polygon by the frozen elevated rims may in this situation slowly flow through the active layer of the rims themselves, from one polygon to another or from polygon center to the ice-wedge channel network (Boike et al., 2008).This phenomenon is important for the water balance: in this case we assume a constant drop of the water table because of lateral runoff R: where R is computed in mm day −1 .Physically it represents both surface and subsurface runoff.We also consider that, due to soil porosity and water holding capacity, the water table rises differently inside and above the soil.We assume drainable porosity s to be 0.7 if S − W t < 0, and 0 if S − W t ≥ 0. We base this assumption on field studies (Helbig et al., 2013;Langer et al., 2011b).
We integrate the model for 90 days, which is approximately the duration of the thaw season in northern Siberia (July-September).when the surface is covered by snow.We also do not consider early summer, when the thaw depth is still limited.External climatic forcings are precipitation and evapotranspiration.We can vary precipitation intensity and its seasonal amount by tuning parameter R p (Eq. 9).In our simulation, we use R p = 30 mm day −1 for a standard scenario, R p = 10 mm day −1 for a dry scenario, and R p = 60 mm day −1 for a wet scenario.For each scenario, we perform ensemble simulations, calculate the average of 30 ensemble members, and compute water table dynamics.We finally derive variations on the amount of landscape area covered by different surface characteristics (i.e., moist centers, saturated centers, or wet centers).
Many studies showed the relationship between water table height and methane emissions, such as Couwenberg and Fritz (2012); Zona et al. (2011); Sachs et al. (2010), andStrack (2004).In order to give an example of the influence of water table and surface wetness on methane fluxes accounting for different climatic forcing, we multiply the area covered by moist, saturated, or wet centers by values that represent the average August methane emissions from each surface type (Sachs et al., 2010).Average methane flux (in mg d −1 ) from the total region at each daily time step is computed following the equation where σ = 94.1 mg d −1 m −2 , ω = 44.9mg d −1 m −2 , and δ = 7.7 mg d −1 m −2 .Values are taken from Table 4 of Sachs et al. (2010).Area C i of polygon centers with different wetness and area r of relatively dry rims are computed in m 2 .In this simple model, we assume polygon rims and moist centers to have the same emission coefficient δ.

PVD and percolation theory
Mathematical properties of PVD are also well suited to describe some landscape-scale physical processes that have been observed in the field.In particular, we focus on interconnectivity properties of the graph, applying percolation theory on PVD.When water flows from polygon centers through the unfrozen top layers of the rims (lateral runoff), it starts to flow into channels on top of the ice wedges between the polygon rims (visible in Fig. 1).Due to thermal erosion, water does not stay confined in the channel but flows through a path of interconnected channels.If this path of channels is large enough to connect polygons from one side of the landscape unit of consideration to the other one, we can assume that water in the channel is able to flow from one side to the other one.The system tested in this study is an island.Consequently, water that was previously confined in the channels would likely run off the island to the Lena river nearby.In our model we assume interpolygonal channels (polygon edges in the PVD, or bonds, in percolation theory) to be active if the polygon experiences lateral subsurface runoff.We also as-sume the polygon channel system to be flat, or with a very small topological gradient, so that outflow due to landscape slope is negligible.
Percolation theory describes the behavior of connected clusters in a random graph: since its introduction in the late 1950s (Broadbent and Hammersley, 1957), it has been widely and successfully applied to a large number of fields, from network theory to biological evolution, and from galactic formation to forest fires (Grimmett, 1999, and references therein).Originally, the question to be answered focused on transport in porous media.Let us imagine a porous stone touching some water on a side: what is the probability that the water passes through the whole stone, coming out at the other side?This simple question led to a basic stochastic model for such a situation and gave birth to percolation theory.In our model, the question would be: what is the probability that an open path of interconnected channels is formed?This is in fact the probability that a giant cluster of channels with flowing water appears.On the PVD, bonds connecting two points of the graph represent channels.Water flowing through rims enters the channel network generated above the ice wedges between the polygons.These channels are underlain by ice-wedges, and water flowing through them may also cause thermo-erosion because of its mechanical and thermal effects on flowing water on ice.In our simple model, if a channel is filled by water, we call it active.Since we assume the environment not to have any slope, all edges of a polygon with W t < TD become active at the same time.
Percolation theory focuses on the probability θ (p) that a giant cluster exists.This probability is nonzero only if the probability p that a bond is active exceeds the critical threshold values p c .This value depends on the graph.Recent studies (Becker and Ziff, 2009) determined numerically the bond percolation threshold p c on a PVD: In order to reach percolation, the condition W t < TD must be reached in a large enough number of polygons.We simulate this result by keeping the thaw depth constant and increasing the water table.We do not consider any physical process involved in water table change, but we keep increasing its level until we reach the percolation threshold.

Results and discussion
In this section we analyze environment statistics by comparing model outputs with available observations, both from field and from aerial photography.We analyze the geometric characteristics of the polygonal tundra and the ones of PVD.We then present and discuss results from 90 days of simulations, and the water table dynamics in different simulated scenarios.Finally, we show results from our bond percolation realization.Our results are qualitatively different from a mean field approximation approach since mean water table variations are computed considering processes taking place in each single polygon.

Polygon statistics
We are interested in the distribution of areas covered by polygon centers C = (1 − q)A pol .As for the polygon area A pol , we compare the histogram of model output for C with a generalized gamma distribution f A and a 2-parameter gamma distribution.Both distributions are accepted by a Kolmogorov-Smirnov test at a 0.05 level of confidence.Parameters of the distribution are taken from literature concerning the polygon area distribution in a PVD (Okabe et al., 2000, and reference therein).Lognormal, exponential, and normal distributions are rejected by the test for both A pol and C.
We then proceed to compare the model results with available data for the ice-wedge polygons (Muster et al., 2012).We tuned the controlling parameter of the point process λ and the random parameter q (Eq.5) with observed mean values of polygon and polygon center areas.Results show that comparisons between model outputs and data are quite robust: not only are the first two statistical moments comparable but they also display a very similar distribution (as confirmed by a χ 2 test at a 0.05 level of confidence).
This result is central in our approach: it shows that PVD description of the landscape statistics is consistent with observations.In particular, we claim that our approach is able to consistently capture the topographic features of the landscape we analyze.Even though the PVD are generated by a completely stochastic process (Poisson point process), they represent area distribution of such a complex environment as the low-center ice-wedge polygonal tundra.In order to test the robustness of this result, we performed simulations on another less random tessellation, i.e., a regular square lat-tice, with polygon area A pol equal to the mean area from observation data (136 m 2 ).Again, polygon center area C = (1 − q)A pol , with parameter q defined in the same way as for the random case.Even though mean polygon and polygon center areas have been tuned with observations, results do not fit polygon center distributions and fail in particular to capture the skewness of the distribution, being therefore rejected by Kolmogorov-Smirnov test.

Response to climatic forcing
Precipitation and evapotranspiration are the main inputs and outputs that drive water level dynamics as described in Eq. 7.
In order to test the model response to drier and wetter conditions, we simulate different configurations by changing the amount of incoming precipitation.We then compare the cumulative precipitation simulated by the model and the one measured in the field in order to have realistic scenarios comparable with data from field campaigns (Boike et al., 2008;Sachs et al., 2008).Field studies observed in particular the influence of lateral runoff and storage: water table drops in the late summer season since the water table usually lies above thaw depth, and it leads to lateral fluxes of water in the ice-wedge channel network.The model captures this process (Fig. 5) as well as the water table variation magnitude.
Model results show that a decrease in precipitation will lead to a drastic drop in the water table level due to lateral runoff (Fig. 5b).Precipitation income is the only input in Eq. 7. Scenarios of such dry seasons are likely to lead to drastically drier soil conditions in the polygonal tundra landscape if there is not a parallel decrease in mean surface temperature, and consequently in thaw depth behavior.On the other hand, an increase in precipitation with respect to the reference value used in the standard scenario will not cause a flooding of the landscape.In the wet scenario we perform simulations with a mean seasonal increase in precipitation of 100 mm yr −1 with respect to the reference value.In agreement with observations, mean water table changes its position by lying slightly above the average polygon surface (i.e., W t < = 10 cm).
Depending on W t position in respect to the surface of the polygon center, soil characteristics in polygon centers vary from moist to saturated and wet conditions.Field studies in this region showed a link between water table position and emissions of greenhouse gases, with respect to methane (Kutzbach et al., 2004;Sachs et al., 2010) and latent heat fluxes (Langer et al., 2011a;Muster et al., 2012).In particular, methane fluxes vary up to an order of magnitude among different surface types.Therefore dynamics of water table also influences dynamics of GHG emissions.Once again, temperature variations and influences on thaw depth seasonal behavior are neglected.If mean surface temperature rises, we expect lateral runoff to start earlier in the season and therefore to be able to drain more water away from the environment.In our wet scenario precipitation further increases.In this case, runoff is not sufficient to balance the excess water, and the area covered by wet centers increases.
Our model is also able to calculate the seasonal dynamics of the landscape area covered by moist centers, saturated centers, and wet centers.These results are dependent on thaw depth dynamics, which we here assume prescribed.Overall, in the ensemble averages of simulations in the three different scenarios we used for precipitation input, the fraction of area covered by saturated centers shows little seasonal dynamics if compared to the fraction of surface covered by wet and moist centers.According to our model, during the summer only extreme drops in water table, or very wet summers, would lead to significant modifications of the fraction of the landscape covered by saturated centers.We argue though, that coupling our simple approach to available permafrost models, and thus with thaw depth behaviors dependent on climatic and environmental conditions, would enable more realistic predictions in different scenarios.In particular, in warmer conditions, soil could thaw deeper and earlier in the season, leading to further drops in the water table level and in the fraction of the landscape covered by different terrains.The simple model here introduced for computing methane emissions is directly related to changes in soil wetness.We show the results of ensemble simulations in different scenarios in the second panel of Fig. 5b.In particular, we observe that a decrease in precipitation (dry scenario) leads to a parallel decrease of methane emissions with respect to the standard case, mainly because of spreading of relatively drier tundra, which presents lower methane fluxes than wet and saturated centers.Water table decreases and more and more centers become drier.According to Eq. ( 6), the area of the lower center decreases at the same time, thus increasing the area covered by rims r, which in our model have the same emission properties of drier centers.If precipitation increases, our model predicts a corresponding increase in methane emission from the landscape.In this extreme case, a more elevated water table would lead to an increased number of wet centers, causing a retreat of the relatively drier tundra.Even though area covered by high-emissive saturated centers decreases, in fact, the decrease of surface covered by relatively drier centers is even larger.This situation is likely to increase methane emissions in respect to the standard scenario.In the wet scenario, however, the decrease in the area covered by moist centers is counterbalanced by the decrease of area covered by saturated centers in favor of wet ones.Therefore, methane emissions, even if, as expected, larger in the wet scenario, do not increase significantly.Our results are only qualitative, since we assume only a fixed average value for the seasonal methane emissions.It is worthy to mention that the order of magnitude of the methane emission is the same as the one found in eddy covariance measurements in the study site (Sachs et al., 2010).This finding has been obtained without any fine tuning and with a very simple methane emission model (Eq.15), but it nevertheless proves the goodness of the model.

Percolation threshold
Figure 7 shows a bond percolation realization.We reach the bond percolation threshold by increasing the water table with a constant thaw depth.More and more polygons experience lateral runoff.Consequently, their channels become active and water flows into the channel network.The number of active channels needed to reach bond percolation is given by the theoretical findings described above, and it increases with a rising water table.We expected the system to reach the percolation threshold if the average water table level lies above the thaw depth, but the correlation between these two phenomena is not trivial.In particular, the system does not reach the percolation threshold at the same time as the average water level is above the average surface.In the real system, if an interpolygonal channel is active, i.e., filled with water, water would not stay confined in it, but would spread to other empty channels because of gravity.This means that, as only single polygons show lateral runoff, this drainage water would be distributed to a certain number of other topographically connected channels, and thus the overall increase of water level within the channels would be dampened.If many polygons experience lateral runoff, the number of empty channels would be drastically reduced.The water of the active channels could then spread in fewer channels.Therefore, the average water level within the channel network would increase.In this case, if a physical threshold somehow analogous to the percolation threshold on the PVD would reached, the water that was confined in the network channel is likely to flow out of the system, namely into a near river or sea.This phenomenon may also have a significant impact on the carbon balance of the system since the outflow of water from the channel is one of the most important phenomena for the exchange of DOC (dissolved organic carbon) and DIC (dissolved inorganic carbon).Improving those F. Cresto Aleina et al.: A stochastic model for the polygonal tundra features could be a future development of our approach, and the focus of further studies, more specifically interested at landscape-scale processes of periglacial environment.Nevertheless, a simple application of mathematical properties of PVD is able to capture landscape-scale behavior of the system.In particular, the number of polygons needed to reach outflow can be easily computed through a bond percolation model.At this stage, water can flow from the system into a near river or to the sea and would be lost.Consequences of strong runoff events are ignored in our model, but field observations suggest that thermoerosion and water runoff could be agents for polygon degradation.We argue that the ability of PVD to simulate a physical process through critical thresholds enhances the power and the theoretical predictive capability of our approach.

Summary and conclusions
Our model provides a simple, stochastic, and consistent approach to upscale the effect of local heterogeneities on GHG emissions at the landscape scale.In particular, we show how applications of Poisson-Voronoi diagrams can statistically represent the main geometric and physical properties of the low-center ice-wedge polygonal tundra.
We describe the landscape using Poisson-Voronoi diagrams, tuning model characteristics with field data.We show that the probability density function of modeled polygon center areas is consistent with the one inferred by observations.Our approach is therefore able to statistically upscale the terrain geometrical characteristics.
Dynamical water table variations in respect to the polygon surface are essential to compute methane fluxes, and we compute such variations during the summer season over the whole region, parameterizing processes within the single polygons.Results of water table dynamics agree with observations.In particular, the stochastic model is able to represent the water balance over the whole landscape.The model relates then the water table position to different soil wetness levels.In particular, we are able to compute for each time step the fraction of the landscape covered by moist tundra (moist centers and polygon rims), saturated centers, and wet centers under different climatic forcing.
Our modeling framework enables us to investigate effects of large-scale interconnectivity among ice-wedge channel network.Applications of percolation thresholds on the PVD suggest possible explanations for water runoff from the landscape.
The stochastic parameterization of our model is far from the precision and accuracy of a mechanistic model.Mechanistic models, however, can reach accuracy only at the local scale, whereas our approach accounts for landscape-scale properties and processes.Our presented approach could be further improved by including parameterization of hydrology, characterization of ice wedges, and phenomena such as thermokarst and polygon degradation.In particular, we argue that linking this approach with available permafrost models would enable even more realistic predictions of water table position and the seasonal development of different surface types.Both variables are fundamental to consistently predict GHG emissions from polygonal tundra.
This model shows a new approach that could be successfully applied to other environments and ecosystems where local processes and microtopography play an important role and a mean field approximation would fail in estimating large-scale features.In order to achieve this goal, we need to gather information on the microtopography of the environment to which the tessellation approach will be applied.For environments such as peatlands, for instance, a characterization of the surface elevation (namely, the amount of surface covered by lawns, hummocks, and hollows) in respect to the water table is essential.In order to properly model other ecosystems, the point process generating the tessellation, which in the present study is a stochastic Poisson point process, could also be substituted by a more regular one.For such a development, it is necessary to adapt the vertical structure and the hydrology parameterization of the model (which now represent polygonal tundra) to the new system under consideration.
Flow network techniques, along with more information on the interpolygonal channel characteristics gathered in the field, could be useful tools to estimate and predict the water flow in the channel network.
Further work will be focused on linking the stochastic model to existing mechanistic surface models, such as land surface schemes (LSS).We think that using the model as an external module for characterizing surface microtopography will improve the variance of the mechanistic model, leading therefore to a qualitative improvement of the mean field approximation for water table level, energy balance, and greenhouse gas emissions.One way to couple these two models is to use an LSS to compute the area covered by such an environment and the climate forcing.The microtopography model driven by the LSS data will then be able to feed back to the LSS the response of local soil features to such forcing for the given area.
Overall, the general agreement between field measurements and model in this study's results suggests that statistical methods and simple parameterizations, if accurately tuned with field data, could be a powerful way to consider spatial-scale interactions in such heterogenous and complex environments.

Fig. 1 .
Fig. 1.Color-infrared aerial picture of polygonal tundra on Samoylov Island, Lena river delta, Russia.Dry rims appear light gray, open water appears black.Dark gray and reddish areas indicate moist to wet areas.

Fig. 3 .
Fig.3.Cross section of a polygon.Surface is divided in drier elevated rims and in lower wet centers.Both surface level and single polygon depth are randomly generated.Water table level W t and thaw depth TD are also randomly initialized, but they dynamically variate with external forcing (W t ) and seasonal cycle (TD).Angle α represents the slope of polygon rims.

Fig. 4 .
Fig. 4. Cumulative probability distribution of polygon center sizes.Observations (blue line), and 100 model realizations (cyan lines) are compared.The figure displays area on the x axis.Mean size of center area is 61.19 m 2 in the observations and 61.36 m 2 in the model on.These values are comparable as confirmed by a Student's t test.Distribution and standard deviation of data and model output are also comparable.Red lines represent the confidence bounds for the data cumulative distribution function at a 0.01 confidence level.

Fig. 5 .
Fig. 5. Ensemble averages of 30 model simulations.The graph in panel (a) displays water table variations over time.Assuming no runoff, and storage term equal to zero, water level should correspond to the blue curve (cumulative precipitation minus cumulative evapotranspiration).Field datasets clearly show a drop in water table at the end of the season, due to lateral runoff, when water table lies above the thaw depth.Panel (b) shows water table dynamics in the three simulated scenarios: wet (blue line), dry (red line), and standard (black line).Shaded areas represent the variance of the ensemble simulations, which increases with the increase of the amount of water input.

Fig. 6 .
Fig. 6.Ensemble averages of 30 model simulations in different scenarios.Panel (a) displays seasonal variations of area covered by saturated centers in a wet (blue line), dry (red line), and standard (black line) scenario.Saturated centers represent the major contribution to methane emission in this area, and in our model their total area drops in both dry and wet scenarios.Panel (b) shows modeled total methane emissions from the whole landscape in the three scenarios.Despite a decrease in the number of high-emissive saturated centers, our simple model shows increased methane emission in the wet scenario because of a drop in the area covered by the relatively drier tundra (moist centers and elevated rims).

Fig. 7 .
Fig. 7. Bond percolation realization: increase of W t level and decrease of thaw depth cause interpolygonal flow of water.Water flows into the ice-wedge channel network, generated by crack processes.In our model, channels become active, and polygons with W t > TD are colored in blue.In panels (a) to (d), the thaw depth becomes deeper and deeper.This process happens seasonally, and with the deepening of the thaw depth, more and more channels become active, until a certain threshold is reached (d), when an open path of interconnected active channels appears.The open path is colored in yellow, and it crosses the whole region: water can flow from channels in the middle of the region to its borders as outflow, and it is not kept in the environment any longer.

Cresto Aleina et al.: A stochastic model for the polygonal tundra 191
table variations at each time step (i.e., each day).The time period we are looking at eff , evapotranspiration ET, and lateral runoff R. Parameter s represents the drainable porosity.We compute water table variations for each polygon at each time step.Water table position is then related to surface characterization, as described in the above paragraph.Since we only Earth Syst.Dynam., 4, 187-198, 2013 www.earth-syst-dynam.net/4/187/2013/F.