DYPTOP: a cost-efficient TOPMODEL implementation to simulate sub-grid spatio-temporal dynamics of global wetlands and peatlands

Abstract. Simulating the spatio-temporal dynamics of inundation is key to understanding the role of wetlands under past and future climate change. Earlier modelling studies have mostly relied on fixed prescribed peatland maps and inundation time series of limited temporal coverage. Here, we describe and assess the the Dynamical Peatland Model Based on TOPMODEL (DYPTOP), which predicts the extent of inundation based on a computationally efficient TOPMODEL implementation. This approach rests on an empirical, grid-cell-specific relationship between the mean soil water balance and the flooded area. DYPTOP combines the simulated inundation extent and its temporal persistency with criteria for the ecosystem water balance and the modelled peatland-specific soil carbon balance to predict the global distribution of peatlands. We apply DYPTOP in combination with the LPX-Bern DGVM and benchmark the global-scale distribution, extent, and seasonality of inundation against satellite data. DYPTOP successfully predicts the spatial distribution and extent of wetlands and major boreal and tropical peatland complexes and reveals the governing limitations to peatland occurrence across the globe. Peatlands covering large boreal lowlands are reproduced only when accounting for a positive feedback induced by the enhanced mean soil water holding capacity in peatland-dominated regions. DYPTOP is designed to minimize input data requirements, optimizes computational efficiency and allows for a modular adoption in Earth system models.


Introduction
Changes in the extent of wetlands affect the climate system biogeophysically and biogeochemically. The surface-toatmosphere exchange of energy and water is fundamentally altered over flooded areas (Gedney and Cox, 2003;Krinner, 2003;Moffett et al., 2010) and wetland ecosystems play a disproportionately important role for the atmospheric methane (CH 4 ) and carbon (C) budgets (Tarnocai et al., 2009;Yu et al., 2010). Today, 175-217 Tg CH 4 , or 20-40 % of total annual emissions originates from wetlands (Kirschke et al., 2013) and the spatio-temporal variability of wetland extent exerts direct control on the seasonal and interannual changes in CH 4 emissions and its atmospheric growth rate (Bloom et al., 2010;Bousquet et al., 2006). Changes in the distribution and productivity of wetlands were most likely a major driver for CH 4 variations during glacial-interglacial cycles and millennial scale climate variability during the last glacial period (Spahni et al., 2005;Schilt et al., 2010).
Wetlands (e.g. marshes, swamps) are ecosystems with their functioning adapted to waterlogged soil conditions. This can be linked to seasonal or permanent inundation where the water table is above surface. Peatlands (e.g. mires, bogs, and fens), are a sub-category of wetlands and are formed when accumulation of organic material exceeds decomposition due to waterlogged, anaerobic soil conditions. Organic peatland soils are characterized by an extremely large porosity where typical values are around 0.8-0.9 m 3 m −3 (Granberg et al., 1999), on the order of 100 % higher than in mineral soils (Cosby et al., 1984). This implies a large soil water storage and retention capacity. Peat- lands not only contribute about 40-50 Tg CH 4 -C to annual CH 4 emissions (Spahni et al., 2011), but also store 500 ± 100 Gt carbon (Gt C) (Yu et al., 2010), which corresponds to about a fifth of the total global terrestrial C storage . In contrast to mineral soils, peatlands continue to accumulate C on millennial timescales owing to the extremely slow decomposition rates and the associated long-lasting legacy effects of climatic shifts that occurred even millennia before today (e.g. the disappearance of the Laurentide ice sheet in the course of the last deglaciation).

Published by Copernicus
Accounting for the pivotal role of wetlands for global greenhouse-gas (GHG) budgets, representations of wetland biogeochemical processes are implemented in land models to hindcast observed past variations and predict future trajectories in atmospheric CH 4 and the terrestrial C balance (Singarayer et al., 2011;Spahni et al., 2011;Kleinen et al., 2012;Melton et al., 2013;Zürcher et al., 2013). Dynamic global vegetation models (DGVMs) and terrestrial biosphere models (TBMs), often applied as modules to represent land processes in earth system models, resolve relevant processes to simulate terrestrial GHG emissions and uptake in response to variations in climate and CO 2 (McGuire et al., 2001), while land surface models (LSMs) are applied to simulate biogeophysical processes associated with the interaction between the land surface and the atmosphere. DGVMs, TBMs, and LSMs, hereafter referred to as land models, often rely on a fixed prescribed extent of wetlands and peatlands. However, predictive model capabilities with respect to the spatial distribution of wetlands and peatlands are crucial when applying models to boundary conditions beyond the present-day state, i.e. when their spatial distribution is substantially different from present-day observational data. Also on shorter timescales, the seasonal and inter-annual variability of wetland extent appears key to explaining the observed variations in CH 4 growth rates (Bloom et al., 2010;Bousquet et al., 2006;Kirschke et al., 2013). In other words, predictions of wetland GHG emissions rely not only on the evolution in the flux rates per unit area, but importantly also on changes in the areal extent of wetlands.
The challenge for global model applications with relatively coarse model grid cells is that even the large-scale hydrological characteristics are determined by the unresolved sub-grid scale topography. Diverse wetland extents simulated by current state-of-the-art land models, applied for bottomup CH 4 emissions estimates, underline this standing issue Wania et al., 2013). Different recent efforts to include dynamical wetland schemes into land models (Gedney and Cox, 2003;Ringeval et al., 2012) are founded on the concepts of TOPMODEL (Beven and Kirkby, 1979). This approach was initially developed to dynamically simulate contributing areas for runoff generation in hydrological catchments. It relies on a topographic index representing the "floodability" of an areal unit within a given river catchment. Using this sub-grid scale topography information, TOPMODEL accounts implicitly for the redistribution of soil water along topographical gradients within a river catchment and predicts the area at maximum soil water content. Neglecting the temporal dynamics of water mass redistribution effects through a channel network topology (river routing), the area at maximum soil water content is used as a surrogate for the inundated area fraction, thereafter referred to as f . TOPMODEL-based implementations have proved successful at capturing the broad geographic distribution of wetlands and their seasonal variability (Gedney and Cox, 2003;Ringeval et al., 2012).
Recently, Kleinen et al. (2012) combined TOPMODEL with a model for peatland C dynamics to predict the boreal peatland distribution and simulate their C accumulation over the past 8000 yr (8 kyr). The rationale for their modelling approach is that conditions for peatland establishment and growth are limited to areas where water-logged soil conditions are sufficiently frequent. Peatland distribution is thus co-limited by f , which is simulated by TOPMODEL.
Here, we present the Dynamical Peatland Model Based on TOPMODEL (DYPTOP). It makes use of the TOPMODEL approach to establish a relationship between the water table depth and the flooded grid-cell area fraction. Once established, this grid-cell-specific relationship is represented by a single analytical function and a set of four grid-cell-specific parameters (provided in the Supplement). This function is used to dynamically predict the inundated area fraction f in combination with the water table depth as simulated by a land model. This simplification reduces required input data, enhances numerical efficiency, and facilitates the adoption of dynamical inundation prediction schemes into land models.
DYPTOP combines this inundation model with a model determining suitability for peatland growth conditions to simulate their spatial distribution and temporal change. This is founded on the approach of Kleinen et al. (2012) but includes a set of modifications to resolve the challenge of predicting the observed spatial heterogeneity of the global peatland distribution across the boreal region. In particular, peatland distribution is considered to be limited by the persistency of inundation, rather than its mean. Furthermore, DYPTOP accounts for the feedback between inundation dynamics, peatland establishment, and the modification of the regional hydrology by the distinct hydraulic properties of organic peatland soils. The present model is designed to account for the temporal inertia of lateral peatland expansion, enabling future investigations of the dynamics of peatland shifts over palaeo-timescales and under future climate change scenarios. In addition, the present study extends the scope of Kleinen et al. (2012) to the global scale, attempts to predict the occurrence of peatland soils also in tropical and sub-tropical ecosystems, and relies on plant physiology parameterizations of peatland-specific plants.
DYPTOP is applied here in combination with the LPX-Bern Version 1.2 Global Dynamic Vegetation Model (see Sect. 2). We start with describing the LPX-Bern model structure in Sect. 2, followed by a detailed description of the DYP-TOP model formulation in Sects. 3 and 4, and a description of the experimental setup in Sect. 5. The model code and required input data are provided in the Supplement and on github (https://github.com/stineb/dyptop). In Sect. 6, we demonstrate that this model framework is successful at reproducing key spatial and temporal characteristics of the dynamics of inundation areas and peatlands on the global scale. These results are discussed in Sect. 7.

The LPX-Bern DGVM
DGVMs simulate processes of vegetation dynamics and terrestrial biogeochemistry in response to climate and atmospheric CO 2 and account for the coupling of the carbon (C) and water cycles through photosynthesis and evapotranspiration. Plant functional types (PFTs) are the basic biological unit and represent different life forms (grasses, trees, mosses) and combination of plant traits (needle-leaved, broad-leaved, etc.). The distribution of PFTs is simulated based on a set of bioclimatic limits and by plant-specific parameters that govern the competition for resources. Here, we apply the LPX-Bern version 1.2, a further development of the LPJ-DGVM (Sitch et al., 2003). It accounts for the coupled cycling of C and nitrogen (N), whereby net primary productivity (NPP) is limited by the availability of explicitly simulated inorganic N species following Xu-Ri and Prentice (2008).
Each grid cell in LPX-Bern is separated into fractions representing different land classes (tiles) where C, N, and water pools are treated separately. Upon any change in the tiles' fractional area, water, C, and N are reallocated, conserving the respective total mass (see Strassmann et al., 2008;Stocker et al., 2014). Here, we explicitly distinguish between natural land with mineral soils (f mineral ), peatlands (f peat ), and former peatlands, now treated as mineral soils (f oldpeat ). We describe how the simulated soil water balance and a diagnosed inundation area (see Sect. 3) can be used to predict changes in the fractional areas f mineral , f peat , and f oldpeat in response to changes in climate and CO 2 (see Sect. 4). Depending on the model application, LPX-Bern land classes may additionally distinguish between land with primary vs. secondary vegetation, croplands, pastures, and built-up areas (see Stocker et al., 2014). These model features are not activated in this study. Ringeval et al. (2014) applied an alternative version of LPX-Bern (version 1.1) to simulate separate C dynamics on floodplains which are represented by a separate land class (tile). This feature is not used for the present study as the focus here is on the spatial dynamics of peatlands and any additional grid-cell tile comes at a substantial computational cost.
Biogeochemical processes and the water balance are simulated using distinct parameterizations on the different gridcell tiles f mineral and f peat . All parameterizations and parameters are identical for f mineral and f oldpeat . On f peat , we apply a version of the LPJ-WHyMe model (Wania et al., 2009b), adopted and modified as described in Spahni et al. (2013). This model simulates peatland-specific soil carbon dynamics that are governed by variations of the water table position and soil temperature. Peatland vegetation is represented by Sphagnum (moss) and Graminoids (sedges). Key parameters such as the decomposition rate of soil organic matter are tuned by Spahni et al. (2013) to best match observational site data (Yu et al., 2010) for peat C accumulation rates over the past 16 kyr. These parameter values are left unchanged for the present study. In contrast to earlier studies of Spahni et al. (2011Spahni et al. ( , 2013, we include three additional PFTs on peatlands. These inherit properties of the tropical evergreen and tropical deciduous tree PFTs and the C4 grass PFT (see Sitch et al., 2003), but are adapted for flood tolerance (Ringeval et al., 2014). Additionally, we removed the upper temperature limitation of the other peatland-specific PFTs, already used in previous studies (Graminoids, Sphagnum) to permit their growth outside the boreal region. Representations for the interaction of the C and N cycles are implemented in the peatland-specific model part as described in Spahni et al. (2013). However, we updated the prescribed soil C : N ratio for the peatland PFTs (C : N ratios for sedges = 35.87, sphagnum moss = 82.35, other woody PFTs = 50.98) according to Loisel et al. (2014).
Parameterizations and parameter values applied for C and N cycling on natural land on mineral soils (f mineral and f oldpeat ) are largely identical to previous applications of LPX-Bern version 1.0 Spahni et al., 2013). Changes since version 1.0 include the application of an improved litter decomposition parameterization following Brovkin et al. (2012). Additionally, the temperature governing soil organic matter decomposition in LPX-Bern version 1.2 is computed based on the simulated temperature profile (instead of a single value representing 25 cm depth, Sitch et al., 2003), weighted by a logarithmic soil C profile, fitted to decreasing C density with depth as measured by Wang et al. (2010) on forest, grass, shrub, and desert ecosystems. Steps 1-3 determine the inundated area fraction f and are described in Sect. 3. Steps 4-6 determine the peatland area fraction f peat and are described in Sect. 4.

Topography and inundated area fraction
TOPMODEL (Beven and Kirkby, 1979) makes use of subgrid-cell scale topography information to relate the grid-cell mean water table position (or water deficit as formulated in the original paper) to the area fraction at soil water saturation within each grid cell. The basic information to determine this  Figure 1. Overview of information flow of DYPTOP. Boxes represent spatial fields of the respective variables, given at the spatial resolution as indicated in the lower right edge of each box. (1) Compound topographic index (CTI) values are derived from the ETOPO1 (2013) high resolution topography data set using the R library "topmodel" (Buytaert, 2011). (2) Fit parameters (v, k, q) are derived by applying the leastsquares fitting algorithm "nls" in R (R Core Team, 2012) to best reproduce the "empirical" relationship between the water table position ( ) and the flooded area fraction (f ) as derived from the ETOPO1 data and Eqs. (2) and (3). f max is the maximum area fraction that is allowed to be flooded within a grid cell and is computed by using a globally uniform threshold value for CTI (CTI min ) below which flooding is prohibited (Eq. 4). (3) (v, k, q, f max ) are prescribed to LPX-Bern to predict f as a function of , which is interactively simulated in LPX-Bern. (4) The potential, hydrologically viable, peatland extent f pot peat is set to the minimum of the N months with highest inundation over the preceding 31 years (see Eq. 13). (5) Peatland C balance criteria and the ratio of precipitation over actual evapotranspiration (POAET) are used to determine whether a peatland can establish in the respective grid cell. (6) If criteria are satisfied, the actual peatland area f peat fraction converges over time to f pot peat . (7) The presence of peatlands affects the local soil water storage and retention capacity and thus exerts a feedback via the mean grid-cell water table position and f . relationship is provided by the sub-grid scale distribution of the CTI. In the following, we refer to "pixels" (index i, here ∼ 1 km) as the grid cells within each model grid cell (index x, here 1 • × 1 • ). The CTI determines how likely a pixel is to get flooded ("floodability"). The higher the value, the higher its floodability. It is defined as where a i represents the catchment area per pixel i, i.e. the total area that drains into/through the respective pixel. β i refers to the slope of the pixel. Here, CTI values are derived from the ETOPO1 high resolution (1 arc min) topography data set (ETOPO1, 2013) and are calculated using the R library "topmodel" (Buytaert, 2011) (Step 1 in Fig. 1). Deriving CTI fields from a topography data set instead of relying on available CTI products allows us to extend CTI fields to areas below the present-day sea level for applications on palaeotimescales.
Following the TOPMODEL approach, we calculate the threshold CTI * x in each grid cell x, as a function of the gridcell-mean water table position x . Here, x is in units of mm above the surface. All pixels i with CTI i > CTI * x are at maximum soil water content. Here, this is interpreted as the respective pixel being flooded. CTI * x is defined by where CTI b is the arithmetic mean CTI value, averaged over the entire primary catchment area b in which the respective pixel is located. This is a simplification in case two pixels i and j exist where CTI i > CTI j , and i lies upstream from j . In this case, the relative floodability of CTI i is affected by the fact that CTI j has a low floodability (low CTI value), when in effect there is hardly any influence as CTI j lies downstream from CTI i . However, CTI values generally increase downstream (drainage area a increases), hence CTI i > CTI * > CTI j is not frequent. Note that the catchment area may extend beyond the model grid cell in which the pixel is located. The catchment area data set is from HYDRO1K (2013). Thus, whether a pixel is flooded, hence CTI i > CTI * x , depends on the floodability of other pixels in the same catchment area. M is handled here as a free (and tunable, see Table 1 and Sect. 7.1.1) parameter. More strictly, M describes the exponential decrease in soil water transmissivity with depth (see Beven and Kirkby, 1979).
Accounting for the full topographical information contained in the CTI values within a grid cell x, the flooded area fractionf x within the respective grid cell is determined by the total area of all pixels within grid cell x with CTI i being larger than CTI * x and larger than CTI min . CTI min represents a lower threshold for flooding, irrespective of the water table and hence for the maximum inundated area fraction in grid cell x: A x is the total surface area of grid cell x, and i runs over all pixels located within the respective grid cell. The choice of CTI min affects the maximum possible extent of inundated land within a grid cell and is further discussed in Sect. 7.1.1.
The distribution of CTI values within a given grid cell and the catchment mean CTI determines the inundated area frac- tionf x of this grid cell for each x (see Eq. 2). This relationship is distinct for each grid cell and is illustrated in Fig. 2 for two example grid cells. Having to rely on the full information CTI i is computationally costly due to the (necessarily) high spatial resolution of CTI i .

B. D. Stocker et al.: DYPTOP
Instead of fitting a gamma function to the distribution of CTI i , as has been done earlier (Sivapalan et al., 1987), we directly define a function fitted to the "empirical" rela-tionshipˆ betweenf and .ˆ is established by evaluatinĝ f using Eqs. (2) and (3) for a sequence of spanning a plausible range of values (here from −2000 mm to 1000 mm) and for each grid cell x individually.ˆ generally has a shape as displayed for an example grid cell in Fig. 2 (black curve) and can be approximated by an asymmetric sigmoid function (red curve) with three parameters (v, k, q) x and for CTI min = 0. Here, we apply monthly mean values of as computed by LPX for each month m and grid cell x: We determine parameters (v, k, q) x using the non-linear least-squares fitting algorithm "nls" in R (R Core Team, 2012) (Step 2 in Fig. 1). Note that x is distinct for each grid cell x, as reflected by the values of parameters (v, k, q) x . Additional information is contained in the cut-off value CTI min that determines the maximum flooded area fraction f x max (see Eq. 4). It follows (Step 3 in Fig. 1) that: The two-dimensional (longitude × latitude) fields (v, k, q, f max ) x carry the full information, here on a 1 • × 1 • resolution, and can be used as input for any global model to predict the (monthly) value of the area fraction that is flooded (f x,m ) from the (monthly) water table position x,m in grid cell x. These data are provided in the Supplement and may be applied in combination with an implementation of Eqs. (6) and (5). An example code programmed as a subroutine in FORTRAN is also provided in the Supplement.
In conclusion, the concept presented here can be described by the re-mapping of toˆ , where the information The choice of CTI min and M determines the parameter set (v, k, q, f max x ) x . The large array CTI is reduced to optimize computational costs. In Sect. 7.1.1, we describe how M and CTI min are constrained using observation-based data.
The water table position ( ) is the ultimate predictor variable for f and is calculated online by LPX-Bern. How is defined exactly may depend on the nature of the soil water model implemented in the respective DGVM, and results for f thus depend on the DGVM-specific predictions of . The following subsection describes the definition of in LPX-Bern. All results shown in Sect. 6 are to be interpreted with respect to this DGVM-specific definition of .

Definition of the water table position
The grid-cell-mean water table position is calculated as a grid-cell fraction-weighted mean, where f mineral , f peat , and f oldpeat are the grid-cell area fractions as described in Sect. 2.
The simulated inundated area fraction f is governed by model predictions of the water table position . In the model for peatland-specific biogeochemistry, peat is the key variable determining soil oxygen status and organic matter decomposition. It is explicitly simulated as described in Wania et al. (2009b) (their Eq. 22). The definition of peat accounts for the particular hydraulic properties of peatland soils. This tends to constrain seasonal peat variations to generally higher values than mineral through mechanisms of enhanced soil water storage and retention and reduced runoff.
On non-peatland soils (f mineral and f oldpeat ), no explicit variable representing the water table position mineral and oldpeat is available in LPX-Bern. In the following, we define mineral and oldpeat as an index that is suitable for the present application.
On non-peatland soils, the water balance, surface and drainage runoff are modelled by a relatively simple "twobucket" approach based on the original LPJ (Sitch et al., 2003). The change in water content of the upper layer is given by the balance between precipitation, snow melt, runoff, percolation to the lower layer, evaporation, and the fraction of plant transpiration extracted from this layer. The change in the lower layer results from percolation from the upper layer, losses to ground water, and transpiration. The soil water model version used here has been extended to account for heat diffusion, melting and thawing across eight soil layers, while the soil water content in the two buckets is uniformly distributed within the upper and lower four layers, respectively. Soil moisture -the governing variable for plant water status -is simulated as a scalar index for each bucket (see Eq. 9) as described in Sitch et al. (2003). This "mixed" approach allows for simulating the restriction of percolation when frozen soil layers are present while still maintaining computational efficiency (as compared to a model where the full water budget and its vertical diffusion/percolation are resolved for each soil layer).
The soil moisture scalar θ i in bucket i varies between 0 at permanent wilting point W PWP and 1 at field capacity W FC and is defined as In the "two-bucket" model, water in excess of W FC is diverted to surface or drainage runoff. This prevents water storage from exceeding W FC and the pore volume to be fully saturated. Hence, the water table position is limited to remain below a certain level * = W FC /φ z, determined by the soil depth z, the porosity φ, and W FC . This may hinder an application of such models in combination with TOPMODEL, as argued in Ringeval et al. (2012).
Here, the monthly updated "water table position" in mineral soils, mineral,m , is defined as an index consisting of the combination of monthly mean water-filled pore space (W l · z l /φ), the monthly total runoff, and the soil depth, modified by the presence of frozen soil layers: Subscripts m, d, and l represent months, days, and soil layers, N m is the number of days for month m, l 0 is the number of the layer just above the first frozen soil layer counted from the top (surface, l = 1), W l,d is the daily updated soil liquid water plus ice fraction in layer l, z l is the soil layer thickness, φ is the porosity (uniform over depth), and runoff m is the sum of monthly total surface and drainage runoff. z l 0 ,d is the soil depth, reduced to the depth at the upper boundary of the uppermost frozen soil layer, if any is present. Otherwise, it is set to the nominal soil column depth z max = −2000 mm. This mimics the amplified susceptibility to flooding on (partially) frozen soils.
However, Eq. (10) may overestimate flooding when the liquid soil water above the uppermost frozen soil layer l 0 is low. Therefore, we applied an effective soil depth z * l 0 ,d instead of z l 0 ,d in Eq. (10), defined as θ d is the daily updated soil moisture index (Eq. 9), averaged over all soil layers above l 0 , and λ is a parameter, here set to 2 (see Sect. 7.1.1). Equation (11) guarantees that at high soil moisture, the effective soil depth z * l 0 ,d is equal to z l 0 ,d . At low soil moisture, z * l 0 ,d is not affected by frozen soil layers and the effective soil depth is always equal to the nominal soil column depth z l 0 ,d = z max = −2000 mm. Lateral expansion and contraction of peatland areas are simulated dynamically as a convolution of (i) peatland carbon (C) balance conditions as simulated by LPX and (ii) flooding persistency as simulated by the TOPMODEL implementation. Peatland C balance conditions are simulated for an area fraction f min peat = 0.001 % in each grid cell globally. This value is small enough not to significantly affect the global C balance (0.0005 % of global peat C according to results presented in Sect. 6), but large enough to provide an effective "seed" for peatland establishment and expansion once conditions for peatland establishment are met (It takes 1158 years from f min peat = 0.001 to 1 at 1 % yr −1 expansion rate, see Sect. 4.3). On this minimum area, we apply the peatland-specific model for C dynamics and the water balance as mentioned in Sect. 2.

Peatland establishment criteria
In each simulation year, a hierarchical series of conditions for peatland expansion or shrinking are evaluated for each grid cell, and the boolean variable pt crit is set (pt crit = "true" if the conditions are met and "false" if they are not; see Fig. 3,and Step 4 in Fig. 1). The primary condition is related to the ecosystem water balance, represented by annual total precipitation divided by (over) annual total actual evapotranspiration (POAET). Global peatland occurrence analyses (Gallego-Sala and Charman et al., 2013) have revealed the limiting role of precipitation over equilibrium evapotranspiration. Here, we apply a threshold of POAET * = 1 to limit peatlands to regions with a positive water balance. Simulated actual evapotranspiration is governed by the water table position and varies between 79.5 and 109.5 % of equilibrium evapotranspiration (EET). This follows from the definition given in Wania et al. (2009a) (their Eq. 23). EET is defined after Prentice et al. (1993) (their Eq. 5).
If this first condition is met, C balance criteria suitable for peatland expansion are satisfied either when peatland soil C accumulates with a multi-decadal average rate of more than 10 g Cm −2 yr −1 , or as long as total soil C exceeds the threshold of 50 kg m −2 . All criteria are computed for each grid cell (note that f peat ≥ f min peat for all grid cells) for the current year by averaging the simulated C balance variables and POAET over the preceding 31 years. This is to reduce interannual variability in pt crit , which is driven by interannual variability in climate (a 31 years time series is repeatedly prescribed during the spin-up; see Sect. 5.2).

Potential peatland area fraction
The potential peatland area fraction f pot peat defines the maximum possible peatland extent within each grid cell. f pot peat is approached during peatland expansion in the case pt crit is "true", taking account of temporal inertia (see Eq. 14). It is determined independently from pt crit and captures both the flooding persistency and the seasonal maximum extent of flooding within the respective grid cell (see Step 5 in Fig. 1). The algorithm applied to determine f pot peat can be described as follows. For each grid cell, all monthly values of the inundated area fraction f m of the preceding 31 years are sorted in descending order. The sorting transforms the vector f to f * .
The potential peatland area fraction f pot peat is then defined as where N is a constrainable parameter. This procedure accounts for inundation persistency as a determining factor for peatland extent, i.e. f * N defines the area fraction that is inundated at least N months during 31 years. We investigated f pot peat , applying values N = (10, 15, 18, 20, 25) (see Sect. 7.2.2). Figure 4 illustrates the sorted vectors f * , for two regions.

Lateral expansion and contraction
Finally, the actual peatland area fraction f peat is simulated as a convolution of pt crit and f pot peat . During the transient simulation (after model spin-up), the annually updated f peat (t) gradually expands towards f pot peat as long as pt crit is "true", and contracts to f min peat , when pt crit is "false". To account for inertia in lateral peatland expansion and contraction, the relative areal change rate is limited to 1 % yr −1 (r = 0.01 yr −1 ).
Upon peatland contraction, the area f peat (t − 1) − f peat (t) is allocated to f oldpeat , and expanding peatlands first expand into f oldpeat . This guarantees that C and N mass on grid-cell area fractions that have never (in the course of the simulation) been covered by peatlands are kept track of separately, and prevents C, N, and soil water from being redistributed across the entire grid cell. At any given time t during the simulation, f oldpeat (t) is thus determined by the maximum peatland area fraction in all preceding years in each grid cell x individually: In LPX-Bern, the monthly varying inundated area fraction is used not only to derive annually varying f peat , but also to simulate monthly varying contributing areas for methane emissions from inundated mineral soils (f inund ). No results for simulated methane emissions are presented in this paper. While contributing areas for methane emissions from peatlands are constant within one year and equal to f peat , f inund is defined by In LPX-Bern, f inund does not affect the C balance on mineral soils, and neither f nor f inund is treated as a separate tile (grid-cell land class).

Peatland water table position feedback
f pot peat is determined by the grid-cell-mean water table position . The water table position on mineral soils is different to that on peatlands for identical driving forces (precipitation, temperature, light) due to different soil properties and different vegetation cover (see Sect. 3.2). Hence, f pot peat depends on the actual peatland area fraction f peat . To illustrate this effect, we additionally plot f * "before" peatland establishment (Fig. 4, left), where f * is determined by the water table position on mineral soils only ( mineral ). This effect is also illustrated by f pot peat,0 vs. f pot peat in Figs. 8 and 9. f peat thus imparts a positive feedback via and the flooded area fraction f through mechanisms of enhanced soil water storage and retention and reduced runoff. Under transiently changing climatic conditions, this leads to a hysteresis behaviour: once peatlands are established, they can persist even under conditions where no new peatlands would form.

Model spin-up procedure for peatland area fraction
Due to the slow turnover times of soil organic matter, pool size equilibration under given environmental conditions is attained only on timescales of 10 3 years for mineral soils and around 10 4 years for peatland soils. Instead of running the model forward over several millennia, we apply an analytical solution to shorten the model spin-up. Equilibrium soil C and N pool sizes (C * ) in models with first-order decay kinetics are defined by their inputs by litter fall (I ), and their turnover times τ : This pool equilibration is applied in spin-up year 1000 for mineral soil pools by averaging I and τ over the preceding 31 years.
Complete equilibration of pools cannot be applied for peatlands due to their turnover times being on the same time scale as their age since initiation. The peatland-specific model spin-up is divided into three phases. Pool sizes are initialized to be empty. In the first phase (here, spin-up years 1-999), the soil and litter C and N pools gradually but slowly increase in response to litter inputs. At the end of phase one, the soil pools are scaled up to near-equilibrium. We assume that present-day litter inputs have been sustained for the preceding t * = 12 000 yr and analytically calculate the respective peatland soil pool sizes as Before this near-equilibration and 200 yr thereafter (second phase), the actual peatland area fraction is held minimal (f peat = f min peat ) in all grid cells. At spin-up year 1200, peatland occurrence conditions (pt crit , see Fig. 3) are assessed in all grid cells and the actual peatland area fraction f peat is directly set to f pot peat , where pt crit is "true" while the temporal inertia of expansion takes no effect. All pool sizes per unit area are held constant at the point of this areal up-scaling and mass is thus not conserved. During the remaining 300 year spin-up time (third phase), temporal inertia and mass conservation are accounted for as during the transient simulation phase. The temporal dynamics of peatland expansion and contraction described in Eq. (14) apply only to the third spin-up phase and the transient period of the simulation, i.e. after the model spin-up.
This spin-up procedure ensures that mineral soils are fully equilibrated, while peatland soils with long turnover times continue to slowly increase in size by the end of the spin-up.

Simulation protocol
Coupled C and N dynamics and the soil heat diffusion and water balance in terrestrial ecosystems are simulated by LPX-Bern, Version 1.2 . This model version is extended to include the DYPTOP model as described in Sects. 3 and 4. Standard parameters are chosen as discussed further in Sects. 7.1.1 and 7.2.2: M = 8 in Eq. (2), CTI min = 12 in Eq. (4), N = 18 in Eq. (13), and λ = 2 in Eq. (11).
Two model simulations were carried out. In the first (S0), peatlands are not accounted for (peatland area fraction is zero everywhere and at all times). In this simulation, the inundation fraction f does not affect the carbon dynamics, nor any other model state variable. In this simulation, = mineral and the potential peatland area fraction before peatland establishment f pot peat can be quantified. In the second simulation (S1), peatlands are accounted for and f is used to determine the peatland area fraction following the method outlined in Sect. 4. In this simulation, is calculated as the grid-cell area fraction-weighted average in mineral and peatland soils (see Eq. 8), and the potential peatland area fraction after peatland establishment f pot peat can be quantified.
For the simulation with peatlands, we apply a spin-up as described in Eq. (18). During spin-up, the model is forced by repeated observational 1901-1931 climate from the CRU TS 3.21 data set (Harris et al., 2013), a constant atmospheric CO 2 concentration of 296 ppm (year 1901value, MacFarling Meure et al., 2006, and nitrogen deposition from Lamarque et al. (2011) fixed at year 1901. The transient simulation period covers years 1901-2012 with observational climate, CO 2 , and N deposition from the same sources. Due to the slow response time scales of peatland area and C pools (centuries to millennia) and the rapid climate and CO 2 changes that occurred during the second half of the 20th century, a spin-up under present-day conditions appears less appropriate. Prigent et al. (2007) combined satellite data from passive microwave, active microwave (scatterometer), altimetry, and Advanced Very High Resolution Radiometry (AVHRR) into a "multisatellite" method to estimate monthly inundated areas over multiple years and covering the entire globe. The updated data set by Papa et al. (2010) is applied here and covers years 1993-2004. This is the first and -to date -only data set that represents the seasonal and inter-annual dynamics of inundation areas. The original data are on a 0.25 • × 0.25 • spatial resolution (at the Equator) and have been regridded for the present application using area-weighted averages (see Fig. 5). Hereafter, "GIEMS" refers to the data set by Papa et al. (2010), which is based on Prigent et al. (2007).

Inundation area
This data set provides information on the temporal variability of inundation that compares well with related hydrological variables (Prigent et al., 2007). However, compared with static wetland maps, the satellite-derived data set of GIEMS notoriously underestimates the inundated area fraction in regions with small and dispersed flooding that amounts to less than about 10 % of the grid-cell area (Prigent et al., 2007). A comparison of GIEMS inundation areas with the Global Lakes and Wetlands Database (GLWD, Lehner and Döll, 2004) suggests that areas classified in GLWD as peatlands ("Bog, Fen, Mire"), "wetlands", and "Swamp Forest, Flooded Forest" are generally under-represented in GIEMS. This mostly affects regions in boreal Canada, Eastern Siberia, Western Amazonia, Congo, and the Tibetan Plateau. This is confirmed by a study focusing on the Amazon catchment and relying on synthetic aperture radar in combination with airborne videography (Melack and Hess, 2010). This regional data product suggests higher inundation area fractions than other remotely sensed data (∼ 15 % averaged over the Amazon catchment). Detecting surface water under dense vegetation generally appears to be challenging due to microwave signal attenuation. Tarnocai et al. (2009) mapped soils in permafrost regions across the northern circumpolar region. For the present study, we converted this data set to a gridded field so that the fraction within each 1 • × 1 • grid cell covered by either histels (peatland soils in permafrost regions) or histosols (peatland soils in non-permafrost regions) defines the distribution of the peatland area fraction. Note that the categorization applied by Tarnocai et al. (2009) reflects the predominant soil type within a given polygon and cannot be directly interpreted in terms of fractional area within a grid cell covered by this type. However, as these data resolve spatial pat-  Fig. 7) is subtracted from the simulated inundation area fraction for a better comparison with GIEMS. Dashed red lines represent simulated inundation additionally corrected for snow cover (areas with with snow cover depth > 30 mm snow water equivalents are masked out) and rice cultivation areas (using the maximum of simulation inundation and wet rice cultivation area after Leff et al. (2004) with f = max(f rice , f )).

Peatland distribution
terns at a high resolution (relying on maps of 1 : 250 000 to 1 : 3 000 000 scale), this transformation appears pragmatic. The same issue applies to the alternative peatland distribution benchmark data set by Yu et al. (2010). These authors provide a map that delineates "peatland-abundant" regions, i.e. where peatlands cover at least 5 % of the landmass. Original binary data on 0.05 • × 0.05 • are regridded here to represent the fractional area with significant peatland cover fraction on the 1 • × 1 • grid applied for the present simulations. This information is not directly comparable to the fractional peatland area but should help here to visualize the global distribution of peatland-dominated regions also in areas outside regions affected by permafrost. The peatland area fraction benchmark data sets are illustrated in Figs. 7, 8, and 9.

Inundation areas
Simulation results suggest that major seasonally inundated areas can be found at high northern latitudes in the Canadian and Siberian tundra with values of f around 25 % and along major rivers in tropical and sub-tropical regions (western Amazon, Ganges/Brahmaputra, Fig. 5). The location and extent of these major simulated inundated areas agree well with observational data (GIEMS), but are underestimated in regions where wet rice cultivation is abundant as rice culti-vation is not accounted for in the present simulations (south and east Asia).
On peatlands, the water table is generally below the surface, which implies that remotely sensed data do not detect or underestimate inundation areas in regions dominated by peatlands. Indeed, the GIEMS data set suggests no significant inundation in regions dominated by peatlands.
Wetland fractions f of around 10 % are simulated in areas of eastern Siberia, the Tibetan Plateau and across large areas of the Amazon basin. These extensive areas of seasonal inundation are not seen in the GIEMS data set. More spatially confined wetland areas with high seasonal maximum values of f across the South American and African continents are captured by DYPTOP, although simulated fractions are lower than suggested by the GIEMS data. Simulated extensive inundation areas in forested regions of the Amazon and the Siberian boreal zone are not captured in the GIEMS data set, while high values in the GIEMS data along water bodies (e.g. Amazon) are not simulated by DYPTOP. Figure 5 (bottom) displays the spatial distribution of the observed and simulated month with maximum inundation over a mean annual cycle. This reveals the large-scale patterns of the seasonal inundation regime. In the tropics, inundation seasonality is driven by seasonality in precipitation and thus ultimately by the zonal shift of maximum insolation over the course of a year. This induces the clear zonal pat-

B. D. Stocker et al.: DYPTOP
terning of maximum inundation between the Northern and Southern Hemispheres, a feature well captured by the model.
In the boreal region, inundation seasonality is dominated by the timing of snow melt. The timing of the seasonal maximum is generally simulated too early compared to observational data. This mismatch is most pronounced in North America. A more detailed regional analysis is conducted below.
Most important wetland CH 4 source regions are located in the tropics (Bousquet et al., 2006) and -to a lesser degree -in peatland-dominated areas of the boreal zone. To assess the simulated inundation seasonality in more detail, we thus focus on a set of regions as indicated by the boxes in Fig. 5 (bottom). The spatial domains are selected to group areas characterized by a similar seasonal inundation regime. Figure 6 reveals that the seasonality of inundation, as well as absolute total inundated area over the course of the season, is well captured by the model. In general, the observed seasonal maxima and minima are closely matched. Mismatches in timing are biggest for the seasonal maximum in high northern latitudes (too early maximum extent in NA and SI) and to seasonal minima in tropical regions of the African (AF) and South American (SA) continents, where the simulated rate of inundation retreat after the seasonal maximum is too rapid.
Across regions, there is no consistency as to whether the model overestimates or underestimates total inundated area and differences are likely linked to regional characteristics. For example, in the region comprising India, China and parts of South-East Asia (IC), the model considerably underestimates inundated area, particularly at its seasonal peak. This has to be interpreted with regard to the fact that anthropogenic modifications of the land surface in areas of wet rice cultivation increase the flooded area beyond naturally inundated regions (e.g. rice paddies constructed on slopes). This anthropogenic extension of flooded areas is most relevant in the wet season, while in the dry season, rice paddies are commonly drained, resulting in an amplification of the seasonal amplitude. Additionally accounting for information on rice cultivation areas improves the agreement between modelled and observed inundation areas in region IC (dashed line in Fig. 6).
In boreal regions, simulated inundation is of relatively short duration and occurs during and after the snow melt when soils are still partially frozen and drainage is inhibited. Compared to observational data, the modelled onset and maximum inundation tend to be too early. This mismatch is most pronounced in NA, where also the maximum extent is underestimated. As indicated in Fig. 6 by the blue bars, simulated inundation onset occurs during months where snow cover is still present. The model is formulated so that f may attain non-zero values as soon as the uppermost soil layer is no longer frozen, irrespective of remaining snow cover. In contrast, satellite-derived data of GIEMS suggest no inundation where snow is present by design (Ringeval et al., 2012). This helps to explain the mismatch in simulated and observed high-latitude inundation in early spring (see dashed lines in Fig. 6, regions NA and SI).

Peatland areas
Simulated total peatland area f peat north of 30 • N is 3.2 mio. km 2 . This is somewhat lower than the range of available estimates. Tarnocai et al. (2009) estimated the total peatland area in boreal permafrost regions to 3.6 mio. km 2 . This estimate is lower than the older estimate of 3.88 to 4.09 mio. km 2 by Maltby and Immirzi (1993) and the value of 4.0 mio. km 2 adopted and reported in Yu et al. (2010), both of which include also peatlands in non-permafrost regions. Simulated tropical (30 • S to 30 • N) peatland area amounts to 0.92 mio. km 2 . This is higher than the value of 0.37 mio. km 2 reported in Yu et al. (2010) and 0.44 mio. km 2 reported by Page et al. (2011). Simulated tropical peatland area in South-East Asia is 0.32 mio. km 2 , higher than the estimate of 0.25 mio. km 2 by Page et al. (2011). Southern peatlands (south of 30 • S) are simulated to cover 0.037 mio. km 2 ; less than the value reported in Yu et al. (2010) of 0.045 mio. km 2 .
The global distribution of the simulated peatland area fraction can be compared to the benchmark maps by Tarnocai et al. (2009) andYu et al. (2010) as displayed in Fig. 7. The model successfully predicts the major peatland areas across the globe. According to the benchmark maps, the largest peat complexes can be found in the Hudson Bay Lowland (HBL) and in the West Siberian Lowland (WSL). Both are simulated by the model with area fraction values on the same order as derived from observations. Also smaller spatial features are well captured. The model suggests significant tropical peatland areas in Western Amazonia and on the South-East Asian islands, in good agreement with the map by Yu et al. (2010). However, these authors suggest important peatland areas also in the Tropics and in the Southern Hemisphere (e.g. the Congo Basin, Patagonia), where the model suggests none or only small peatland extent.
In the following, a focus on the two regions where the largest peatland complexes are located will serve to illustrate these model predictions and allow a more detailed comparison with the benchmark maps.
As outlined in Sect. 4, the distribution of the peatland area fraction f peat is simulated as the combination of (i) the suitability of climate and peatland vegetation growth conditions for long-term C accumulation in soils, (ii) the flooding persistency, and (iii) the effect of peatland presence on the regional-scale hydrology by imposing a positive feedback on the extent of peatlands. These three steps are visualized as the potential peatland fraction f Original YU data delineate grid cells with a significant peatland cover fraction (> 5 %). Original binary data on 0.05 • ×0.05 • are regridded to represent the fractional area with significant peatland cover fraction on the 1 • × 1 • grid. This information is not directly comparable to other panels and is therefore illustrated with a distinct colour key. Bottom row: simulated, f pot peat is the potential, hydrologically suitable peatland area fraction after peatland establishment, f peat is the simulated actual peatland area fraction taking account of the carbon balance criteria. tion. Figures 8 and 9 illustrate these three steps for the boreal regions of North America and Siberia.
The spatial distribution of f pot peat reflects both the topography and the soil water balance, suggesting that areas in North America with the highest extent and persistency of flooding are located around the Hudson Bay including large areas on Baffin Island, across a large area of Quebec around 50 • N, and in south-western Alaska. In Siberia, large f pot peat are simulated across the WSL at around 60 • N, the North Siberian Lowland at around 70 • N and 90-110 • E, and along the north-east Siberian coast.
In areas where peatlands are simulated to establish, the mean water table position is generally lifted upwards and flooding persistency tends to be extended. This increases the simulated potential peatland area fraction to values of around 0.9-1.0 along the southern coast of Hudson Bay (HBL) and 0.5-1.0 in the WSL. Outside areas of significant peatland occurrence, this mechanism takes no effect and thus separates peat-dominated areas from their surroundings and results in the high spatial heterogeneity found by Tarnocai et al. (2009). Although peatlands are simulated to establish in the Quebec region of 60 to 80 • W, and 50 • N, f peat does not attain values as high as in the HBL. This is ultimately due to the limit to the maximum inundated area set by the choice of CTI min in Eq. (3). In other words, topographical properties do not allow for extensive peatland establishment as in the flat terrain of the HBL.
Another way to display this effect is visualized in Fig. 4 which illustrates the array of ranked inundation fractions for each grid cell (f in Eq. 12) before (left) and after (right) peatland establishment. In the latter case, inundation is extended throughout the season and affects larger area fractions. Moreover, this mechanism tends to affect mostly those cells that feature large peatland area fractions also according to Tarnocai et al. (2009) and is thus crucial to predict spatially concentrated peatlands in large flatlands.
Other major peatland regions suggested by Yu et al. (2010) around Great Bear Lake (55-55 • N/120 • W) and in Eastern Siberia are under-represented by the model mainly due to topographical restrictions (see f pot peat ). However, benchmark maps are not consistent with respect to the extent and presence of peatlands in Eastern Siberia.
At higher latitudes of the tundra regions, peatland growth conditions (pt crit ) are mainly responsible for limiting their establishment. Model predictions are consistent with the maps of Tarnocai et al. (2009) andYu et al. (2010) in suggesting no significant peatland occurrence beyond a climatical northern frontier where cold temperatures limit plant productivity as illustrated in Fig. 8.
Simulated global scale controls of peatland occurrence are illustrated in Fig. 10. Beyond a southern frontier in Eurasia and the western American continent, peatland establishment is primarily limited by the hydrological balance expressed as POAET. In more humid regions of the temperate zone, as well as tropical and sub-tropical areas, peatland occurrence   Original YU data delineate grid cells with a significant peatland cover fraction (> 5 %). Original binary data on 0.05 • × 0.05 • are regridded to represent the fractional area with significant peatland cover fraction on the 1 • × 1 • grid. This information is not directly comparable to other panels and is therefore illustrated with a distinct colour key. Bottom row: simulated, f pot peat,0 is the potential peatland area fraction, considering hydrological suitability without the peatland-water table position feedback, f pot peat is the potential peatland area fraction, considering hydrological suitability including the peatland-water table position feedback; f peat is the simulated actual peatland area fraction, taking account of the peatland establishment criteria (pt crit ) and peat expansion and contraction. is largely limited by the long-term soil carbon balance. In these regions, the difference between litter inputs (governed by NPP) and decomposition rates (governed by soil temperature and moisture) is not sufficiently large to allow for longterm C accumulation in peatland soils. In the remaining areas, LPX simulates suitable conditions for peatland establishment, but their extent is limited by the topographical setting and ultimately by the simulated inundation persistency. The global overview of Fig. 10 reveals the dominant role of topography to limit peatlands not only along major mountain ranges (e.g. Ural, Rocky Mountains), but also in eastern Siberia and Quebec. Smaller areas with long-term C accumulation in peatland soils are simulated in the mid-latitudes and the tropics, but these appear to be located mainly in areas where topography and inundation persistency limit peatland extent.

Peatland carbon
Simulated global C stored in peatland soils is 555 GtC (mean over years 1982-2012), with 460 GtC stored in northern, 88 GtC in tropical, and 8 GtC in southern peatlands. This is broadly compatible with the estimate by Yu et al. (2010) of 547, 50, and 15 GtC for northern, tropical, and southern peatland C stocks.
Note that C storage in all peatland soils is simulated under the assumption that accumulation occurred over 12 kyr of constant pre-industrial climate and CO 2 (see Sect. 5.2). This simplified setup is chosen to assess the capabilities of a dynamic peatland model without having to rely on information of the climatic past. Therefore, values should not be considered as an explicit estimate for present-day peatland C storage and are thus not highlighted further.

Discussion
The TOPMODEL approach presented here provides a costefficient solution to the challenge of dynamically simulating the global distribution and the seasonal variation of inundated areas. We combine this information with simulated C accumulation in persistently inundated soils to predict the spatial distribution of peatlands and its temporal change.

TOPMODEL implementation
Inundation is constrained to topographically conditioned areas, which must necessarily be treated at the sub-grid scale in any global model. Here, we rely on a TOPMODEL approach to establish a relationship between the soil water balance and the inundated area fraction for each grid cell and describe this relationship using a set of four fitted parameters for each grid cell. These parameter fields are made freely available and can be prescribed to any land surface or vegetation model in com-bination with the dynamically modelled soil water balance to predict inundation extent. This opens up new possibilities to simulate effects of changes in inundation areas on the climate system and enables modelling studies to extend their temporal scope. This is relevant for modelling changes in wetland distribution and associated changes in CO 2 and CH 4 emissions over inter-annual to millennial timescales both for the past and for the future, and to quantify associated climatewetland feedbacks. We tested the model against a remotely sensed data product for the monthly global distribution of inundated areas (Prigent et al., 2007;Papa et al., 2010). However, TOP-MODEL has originally been developed to simulate the area fraction at maximum soil water content (Beven and Kirkby, 1979) and model predictions are therefore not directly comparable to flooding data that represent areas where the water table is above the surface. TOPMODEL predictions for the area fraction at maximum soil water should be regarded as a surrogate for the inundation area fraction that should follow similar spatial and seasonal patterns and exhibit a similar sensitivity to climate change.

Choice of model parameters
Apart from LPX-specific variables related to the soil water balance, the simulated inundated area fraction f is governed by the function and is thus sensitive to the choice of parameters M (in Eq. 2) and CTI min (in Eq. 3). Similarly, the peatland area fraction f peat depends not only on LPX's predictions for the soil C balance, but (in addition to M and CTI min ) also on the choice of N in Eq. (13) and λ in Eq. (11) (for a discussion on peatland-specific parameter choice, see Sect. 7.2.2).
M represents a physically based parameter describing the exponential decrease of soil transmissivity with depth (Beven and Kirkby, 1979). Here, we consider M as a tunable but globally uniform parameter. This is in contrast to Ringeval et al. (2012), who modified the CTI values to obtain best results. We tested the model performance in terms of simulated f for a range of parameter values M and CTI min that yield plausible results for the total simulated inundated area f . Then, given a selected parameter combination (M, CTI min ), we assessed a range of parameter values N to simulate the potential peatland area fraction f pot peat . Increasing M causes a general contraction in f . Note, however, that M and f do not relate linearly, but depend on the distribution of CTI. CTI min "caps" the maximum flooded area fraction in each grid cell and thus limits f in areas with generally low CTI values (mountainous regions). We first constrained CTI min to a range that appears plausible. The effect of varying CTI min within this range (here we assessed CTI min = 10, 12, and 14) is rather small for the annual mean total inundated area but slightly affects the seasonal amplitude.
In a second step, we assessed different parameter combinations (M = (7, 8, 9), with CTI min = (10, 12, 14)) by visually comparing results with observational data from Prigent et al. (2007). Due to the difference in the nature of the observational data set and the model applied here (see also Sect. 7.1), we could not apply quantitative criteria to constrain M, and CTI min . Instead, we relied on a visual comparison and selected a standard choice of M = 8 and CTI min = 12, so that major tropical and sub-tropical wetlands are captured while limiting the overestimation of total inundated area. A documentation of this parameter exploration can be found in Stocker (2013).
In general, none of the parameter combinations resulted in the spatially confined and concentrated spatial pattern of inundated areas suggested by Prigent et al. (2007). An exploration of a broader range of parameter value combinations partly resolved the apparent differences in spatial heterogeneity but resulted in a pronounced underestimation of the seasonal variability and an overestimation of total inundated area.

Comparison with GIEMS
The model is generally successful at capturing the global distribution of the seasonal maximum inundated area fraction and the seasonal timing of maxima across the globe. Differences in observed and simulated maximum inundation are mainly linked to the spatial pattern and the distribution of values for maximum inundation. While the TOPMODEL approach suggests large areas of extensive inundation with relatively low values, the GIEMS data suggest more spatially confined inundation areas and feature areas with high values that are not captured by the TOPMODEL approach. For example, simulated extensive inundation with values around 5-20 % across large areas in the Amazon region, the Tibetan Plateau, or Eastern Siberia appear not to be supported by the GIEMS data. In contrast, high values along rivers (e.g. Amazon, Mississippi, Euphrates, Ganges, Brahmaputra) and in regions containing major or numerous inland water bodies (e.g. boreal Canada, Paraná, Pantanal, Lake Chad) are not captured by our TOPMODEL implementation. This apparent disagreement has to be interpreted with regard to the caveats of the GIEMS data set mentioned in Sect. 5.3.1. A similar mismatch between observations and TOPMODEL-based simulation results was found by Ringeval et al. (2014). Extensive inundation is simulated by DYPTOP in areas classified as "Flooded Forest" or "Wetland" in the Global Lakes and Wetlands Database (Lehner and Döll, 2004). Melack and Hess (2010) quantify the "floodable" area fraction of mapped areas within the Amazon basin at 15 %. This agrees well with the seasonal maximum inundated area fraction across the Amazon catchment of 13 % suggested by our results (average over 1992-2004). Moreover, the model we applied here relies on a land mask that defines the actual land area fraction within each grid cell and thus excludes permanent water bodies. Consequently, the maximum simulated inundation area fraction is limited to the land area fraction in the respective grid cell, while values in the GIEMS data include permanent water bodies. This conceptual difference in the nature of the observational vs. model data contributes to the apparent disagreement in regions with extensive water bodies as mentioned above.
Third, high values of observed annual maximum inundation area in India, around the Mekong river, and in southern and eastern China are due to widespread wet rice agriculture with "anthropogenic wetlands" in the form of rice paddies. The model applied here does not account for any anthropogenic land use or land cover change. It is to be noted that these anthropogenic land surface modifications appear to have resulted in a substantial amplification of naturally occurring flooding and its seasonal amplitude.

Regional characteristics
While hydrological studies commonly focus on the scale of river basins, inundation is not necessarily confined to an individual basin. We thus investigated six deliberately selected continental-scale regions where each region is characterized by a similar seasonal hydrological regime and contains some of the major global wetlands. Within each region, the model broadly captures the observed range of total inundated area and the timing of seasonality. Simulated areas do not exhibit any consistent bias across regions and model-observation differences appear to be linked to land cover characteristics in individual regions (e.g. extensive forest cover, or "anthropogenic wetlands" as mentioned above).
The model applied here neglects the temporal dynamics of downstream water redistribution within the catchment. Instead, inundation is driven by the variations in the immediate soil water balance which does not account for delayed effects of preceding runoff. This aspect is likely to contribute to the overestimated rate of inundation retreat after the seasonal maximum in areas with large river basins (see SA and AF in Fig. 6). Such a flooding hysteresis has also been discovered by Prigent et al. (2007) by comparing precipitation seasonality with inundation seasonality.
In high-latitude regions of North America and Siberia, a similar hysteresis between river discharge and inundation area has been identified by Papa et al. (2008). The model applied here fails to reproduce this pattern, with the seasonal maximum inundation being too early and the retreat too quick. The former may be linked to a crude model representation of snow melt, ice jams in the river valley during early spring (Ringeval et al., 2012), or the fact that inundation is simulated as soon as the uppermost soil layer is no longer frozen, irrespective of remaining snow cover, which would prevent satellites from detecting water. A similar mismatch has been found by Ringeval et al. (2012) in boreal North America. The overestimated rate of retreat may also be related to the structure of different river basins where dis-connections between the river channel and floodplains may cause delayed inundation retreat but it is not captured by the model (Ringeval et al., 2012).

TOPMODEL in combination with soil moisture models
Model predictions for inundation areas are determined by the simulated soil water balance. However, soil moisture across the soil profile, percolation, and runoff generation are often not physically resolved in vegetation and land surface models. This makes it difficult to evaluate modelled soil moisture against observational data, which themselves are subject to notorious caveats (Schumann et al., 2009). In LPX, soil moisture is modelled as an index ranging from 0 at the permanent wilting point to 1 at field capacity, while water in excess of the field capacity is diverted to runoff. This prevents the soil pore volume from being fully water-filled and hence the water table position from reaching the surface. Yet TOP-MODEL essentially relies on the information of the water table depth (or deficit to saturation). How can this challenge be met? Here, we define " " as an index combining soil water content and runoff. This resolves the problem of notoriously low actual water table positions in index-based soil moisture models. Furthermore, is modified to account for the presence of impermeable frozen soil layers. This leads to a higher susceptibility to flooding in affected regions. Additionally, we tested to what degree additional information on the drainage capacity (permeability) of the sub-soil substrate could help to improve simulation results. The new data set on sub-soil permeability by Gleeson et al. (2011) is designed for global applications in land surface/vegetation models. We found that in combination with a soil water balance model of the type implemented in LPX, this additional information does not suit its purpose as soil moisture in the upper layers is hardly affected by drainage from the lowest layer. However, an implementation of this data set may have great potential in combination with a soil water balance model where percolation across soil layers and runoff are simulated based on physical diffusion equations and infiltration limitation (Ekici et al., 2014).

TOPMODEL as a diagnostic
The implementation of the TOPMODEL approach described here can be regarded as a simple diagnostic function of an independent variable, simulated by the vegetation/land surface model (here the index ). This TOPMODEL implementation can also be applied offline as no feedback exists between simulated inundation areas and the soil water balance, runoff, and biogeophysical land surface properties. This is in contrast, for example, to the study by Ringeval et al. (2012) who used the TOPMODEL concept in a global vegetation model to improve runoff predictions. Given that the main focus of the present study is on the TOPMODEL implementation (and how it can be used to predict peatlands) and not on the vegetation model's (here LPX's) prediction of the soil water balance, we do not present any further evaluation of LPXspecific hydrological variables (e.g. soil moisture, runoff). Similarly, we refrain from using simulated CH 4 emissions to constrain inundation area and variability to avoid using indirect information that is subject to its own specific model uncertainties.

Dynamical peatland model
Peatlands may establish where soil conditions are suitable for long-term accumulation of organic matter. Both litterfall and decomposition rates exert direct control on the soil C balance, but the latter may vary most as heterotrophic activity responsible for decomposition is largely inhibited under anoxic conditions, i.e. when soils are waterlogged. However, water saturation/flooding is constrained by the local topographical setting and any prediction of the peatland distribution has to account for this sub-grid scale information. Therefore, we applied a TOPMODEL implementation to account for soil moisture redistribution within a catchment area and to dynamically determine where inundation is sufficiently persistent for peatland establishment, and combined this with a model for C and water dynamics in peatland soils.
Most previous modelling efforts had to rely on externally prescribed maps defining the peatland distribution based on present-day observations, and available paleoecological syntheses have relied exclusively on basal dates of existing peatlands (MacDonald et al., 2006;Yu et al., 2010;Yu, 2011), thus implicitly ruling out the existence of peatlands that have now disappeared. Also the response of peatland extent under future scenarios of climate change has been contradictory (Ise et al., 2008;Gignac et al., 1998;Bragazza et al., 2008), partly owing to the unresolved processes of lateral expansion and retreat. Kleinen et al. (2012) have proposed a solution to some of these challenges by combining a TOPMODEL approach with the LPJ-WhyMe model for peatland C and water dynamics. The DYPTOP model presented here follows the same path and extends the scope by adding the temporal dimension of peatland expansion and retreat in response to changes in climate, CO 2 and (potentially) the presence of ice sheets. This opens up a new perspective on the terrestrial C balance over multi-millennial timescales and glacialinterglacial cycles where peatlands may both appear and disappear in different regions.
Compared to the model presented by Kleinen et al. (2012), DYPTOP also appears successful at producing the previously unresolved spatial heterogeneity of the global-scale peatland distribution. This is mainly achieved by accounting for the elevation of the grid-cell-mean water table depth by peatlands and their large water retention capacity. Additional tests (not shown) have also revealed that it is crucial to average CTI values in Eq. (2) over the respective catchment area, and not just the respective model grid cell as described in Kleinen et al. (2012).

Choice of the simulation setup
For simplicity, the model documentation presented here is focused on a near-equilibrium state of peatland distribution and C storage, representing 12 kyr of sustained soil organic matter accumulation under constant pre-industrial climatic conditions and CO 2 , which mimics the relatively constant pre-industrial Holocene conditions (since 11.7 kyr before present) (Wanner et al., 2008). This simplification limits the direct evaluation to assessing the climate space in which peatlands are simulated to establish and persist today but does not allow for a direct evaluation of the rate of lateral peatland expansion and contraction. In the approach chosen here, peatland area fraction may scale up from the minimum fraction of 0.001 to 100 % in 1158 yr. This is set by the model parameters f min peat , and the relative areal change rate of 1 % yr −1 in Eq. (14). This approach assumes that expansion is proportional to the peatland area and implies exponential areal growth where the potential peatland area fraction is attained on centennial to millennial timescales after initiation (pt crit switched to TRUE). The choice of these parameters does not significantly affect the results presented here as shifts in the spatial peatland distribution were relatively minor throughout the 20th century. Simulated peatland C storage in grid cells not fulfilling establishment criteria (pt crit = FALSE) is only 2.9 TgC (0.0005 % of the global simulated peat C at 1900 (570 PgC)) and is therefore negligible for global C budgets. Further studies could be aimed at assessing these temporal dynamics by benchmarking DYP-TOP driven by transiently changing climate and CO 2 since the Last Glacial Maximum. Peatland initiation could be used as a target variable and compared to observational data on basal ages (MacDonald et al., 2006).

Choice of model parameters
We assessed the simulated peatland area fraction and total C storage for a range of DYPTOP model parameters N (see Eq. 13), λ (see Eq. 11), C * peat , and dC * peat dt (see Fig. 3) and compared results with data from Yu et al. (2010) and Tarnocai et al. (2009). Increasing N reduces f pot peat and vice versa. λ affects mineral and increasing values reduce f pot peat in regions affected by permafrost (most effectively in east Siberian boreal regions). This can be assessed offline, as f pot peat represents the potential peatland area fraction before peatland establishment and depends only on f = ( mineral ). However, the simulated actual peatland area fraction is subject to the effects of peatland establishment and an offline optimization to constrain parameters is not possible. Furthermore, the simulated peatland soil C pools and therefore (indirectly) f peat depend on the full history of C accumulation since peat initiation (∼ 10-15 kyr in reality, MacDonald et al.,  Tarnocai et al. (2009). Original YU data delineate grid cells with a significant peatland cover fraction (> 5 %). Original binary data on 0.05 • × 0.05 • are regridded to represent the fractional area with significant peatland cover fraction on the 1 • × 1 • grid. This information is not directly comparable to other panels and is therefore illustrated with a distinct colour key. Bottom row: simulated, f pot peat,0 is the potential peatland area fraction, considering hydrological suitability without the peatland water table position feedback, f pot peat is the potential peatland area fraction, considering hydrological suitability including the peatland water table position feedback. f peat is the simulated actual peatland area fraction, taking account of the peatland establishment criteria (pt crit ) and peat expansion and contraction.
peatland limitations inundation persistency POAET C balance topography Figure 10. Global distribution of limitations to peatland establishment. The mapping follows from the assessment of pt crit as described in Sect. 4.1. The primary limitation is given by precipitation divided by evapotranspiration (POAET < 1). Long-term C accumulation sets a secondary limitation (purple). Within the green areas, peatland area fraction may extend to its potential maximum f pot peat (not shown), which is limited by inundation persistency. The latter is subject to the soil water balance and topography. Limitation by topography is represented by the theoretical maximum inundation fraction f max (see Eq. 4), shown as red for f max = 0 and green for f max = 1, and interpolated colours for values in between.
2006). A comparison with observational data is thus necessarily confounded by the fact that the present study relies on a schematic 12 kyr spin-up (see Sect. 5.1). Immense computational resources required for transient spin-ups covering the last 10-15 kyr prevented us from conducting a comprehensive parameter exploration. Instead, we applied a coarseresolution setup and tested plausible parameter combinations in a simplified transient spin-up. We selected N = 18 and λ = 2 as a parameter combination that yields a good agreement with respect to the total peatland area (∼ 4 mio. km 2 , Yu et al., 2010) and its distribution across different regions. The fact that this choice simultaneously and reasonably complies with observational estimates in terms of total peatland C mass (365-550 GtC, Tarnocai et al., 2009;Yu et al., 2010) is independent from the parameter tuning deployed here. Instead, parameters and parameterizations of the model for peatland C and water dynamics are left unchanged from an earlier version  where parameter values have been tuned to match site-scale observations of C accumulation since their establishment. However, in order to extend the scope of this earlier study to peatlands outside the boreal region, we introduced three additional PFTs, suitable for growth in warmer climates.
The mass balance criterion dC * peat dt determines whether conditions for long-term peat soil C accumulation are satisfied. This is relevant mostly for peatland initiation (at early stages, the criterion for C * peat is not satisfied). Additional transient long-term spin-ups showed that dC * peat dt = 20 g Cm −2 yr −1 would be too restrictive for North American peatlands to establish (now shown). Our choice of dC * peat dt = 10 g Cm −2 yr −1 is motivated by observational analyses that suggest that the vast majority of examined peats exhibit longterm C accumulation rates above this value (Charman et al., 2013). The C density criterion C * peat is not independent from dC * peat dt as it reflects a time-integration of the latter. That is, after millennia of sustained peat C accumulation, soil properties are sufficiently altered and the land qualifies as a peatland even when dC * peat dt is too low. This is relevant when conditions become unfavourable for new establishment and introduces a hysteresis effect. The choice of C * peat = 50 kg C m −2 is chosen to reflect typically observed peatland soil C contents (Tarnocai et al., 2009). However, the variability is large. Again, the choice of this value is not critical for the results presented here, where the vast majority of peatlands have soil C contents greater than 100 kg C m −2 and no large climate shifts are affecting the peatland distribution.

Global scale controls of peatlands
The modelling of the global occurrence and extent of peatlands presented here relies on relatively simple governing rules and reveals different controls of peatland occurrence across the globe (see Fig. 10). Peatland occurrence is assumed to be primarily constrained to a positive ecosystem water balance, here assessed by the fraction of precipitation over equilibrium evapotranspiration. Within the space set by this limit, C accumulation in peatland-type soils has to be sufficient to build up a significant organic soil horizon.
If these criteria are met, peatlands may establish in the respective area. Peatland extent then depends on the topographical setting, which limits the area that gets flooded sufficiently often to allow for suitable soil moisture conditions inhibiting decomposition of soil organic matter. This approach accounts for flatland-type peatlands but cannot predict other peatland types, e.g. blanket bogs, which appear to be solely limited by an extreme ecosystem water surplus and not by topography (Gallego-Sala and . This is likely to explain the difference in simulated vs. observed peatland occurrence in cool and wet mountainous regions (Pacific coast of Canada and Alaska, Patagonia, United Kingdom and Ireland, Scandinavian Atlantic coast, Alps).
Our analysis also showed the pronounced difference between the potential peatland area fraction before (f pot peat,0 ) and after (f pot peat ) peatland establishment. This is more than a technical aspect but reveals a sponge effect of peatlands on the regional water balance. Their exceptionally large porosity allows more water to be stored in soils, and thus reduces runoff and maintains a higher water table throughout the season. This is crucial particularly in boreal regions, where a size-able fraction of annual soil water input is provided by snow melt. Peatland occurrence thus feeds back to improved conditions for peatland expansion via enhanced water retention. In the model, this feedback arises as the grid-cell-mean water table depth is an area-weighted average of the water table depth in mineral and in peatland soils (red arrow in Fig. 1). This leads to a larger area fraction being flooded throughout the year when peatlands are accounted for (see Fig. 4) and is essential to resolve the observed characteristic spatial pattern of peatland area fractions with sharp contrasts between boreal lowlands (HBL, WSL) and surrounding areas. Future model development may account for altered soil parameters and water retention capacity on f oldpeat due to an elevated soil organic matter content compared to other mineral soils on f mineral . This may add to the hysteresis behaviour of peatlands when conditions become unsuitable for new establishment during transient simulations.

Conclusions
The DYPTOP model presented here incorporates a TOP-MODEL approach using sub-grid topography information to simulate spatio-temporal dynamics of inundation and peatland establishment. The regional total inundation extent and its seasonality agree well with observations, although the stark spatial heterogeneity suggested by remotely sensed inundation data is not fully captured by DYPTOP. This is a common result of TOPMODEL-based models.
DYPTOP successfully predicts peatland distribution across continental scales and rests on an approach where inundation persistency is used as a constraint that is applied in combination with the simulated C balance in organic (peatland-type) soils to simulate the spatial extent of peatlands. We have demonstrated that DYPTOP provides a solution to the challenge of reproducing the characteristic spatial heterogeneity of the peatland distribution by accounting for their sponge-effect on the local water balance. Enhanced water retention allows peatlands to extend and cover large lowland regions.
DYPTOP is designed to minimize input data requirements, optimizes computational efficiency, and allows for a modular adoption of respective code for application in earth system models. This opens up new opportunities to investigate the response of the wetland and peatland distribution, their carbon storage and methane emissions to large climatic shifts as observed in the past and predicted for the future.

Code availability
The code for the DYPTOP model as described in Sects. 3 and 4 is provided in the Supplement and through the openaccess online code repository github (https://github.com/ stineb/dyptop). It is programmed in Fortran 90 and can be compiled using the PGF90 or gfortran compilers (others not tested). This code bundle is designed as a demonstration for how to implement DYPTOP in a global vegetation model and calculates DYPTOP variables for an example grid cell in the Hudson Bay Lowlands (89.5 • W/55.5 • N).
The Supplement related to this article is available online at doi:10.5194/gmd-7-3089-2014-supplement.