The role of hydrological transience in peatland pattern formation

The sloping ﬂanks of peatlands are commonly patterned with non-random, contour-parallel stripes of distinct microhabitats such as hummocks, lawns and hollows. Patterning seems to be governed by feedbacks among peatland hydrological processes, plant micro-succession, plant litter production and peat decomposition. An improved 5 understanding of peatland patterning may provide important insights into broader aspects of the long-term development of peatlands and their likely response to future climate change. We recreated a cellular simulation model from the literature, as well as three subtle variants, to explore the controls over peatland patterning. Our models each consist 10 of three submodels, which simulate: peatland water tables in a gridded landscape; a simple representation of microhabitat dynamics in response to water-table depths; and changes in peat hydraulic properties. We found that the strength and nature of simulated patterning was highly dependent on the degree to which water tables had reached a steady state in response to hydro- 15 logical inputs. Contrary to previous studies, we found that under a true steady state the models predict largely unpatterned landscapes that cycle rapidly between contrasting dry and wet states, dominated by hummocks and hollows, respectively. Realistic patterning only developed when simulated water tables were still transient. Literal interpretation of the degree of hydrological transience required for pattern- 20 ing suggests that the model should be discarded; however, the transient water tables appear to have captured some aspect of real peatland behaviour that generates patterning. Recently-buried peat layers may remain hydrologically active despite no longer reﬂecting

understanding of peatland patterning may provide important insights into broader aspects of the long-term development of peatlands and their likely response to future climate change.
We recreated a cellular simulation model from the literature, as well as three subtle variants, to explore the controls over peatland patterning. Our models each consist 10 of three submodels, which simulate: peatland water tables in a gridded landscape; a simple representation of microhabitat dynamics in response to water-table depths; and changes in peat hydraulic properties.
We found that the strength and nature of simulated patterning was highly dependent on the degree to which water tables had reached a steady state in response to hydro- 15 logical inputs. Contrary to previous studies, we found that under a true steady state the models predict largely unpatterned landscapes that cycle rapidly between contrasting dry and wet states, dominated by hummocks and hollows, respectively. Realistic patterning only developed when simulated water tables were still transient.
Literal interpretation of the degree of hydrological transience required for pattern-1 Introduction

Background
The surface of northern peatlands often comprises a patchwork of distinct, small-scale 5 (< 10 m; known as scale-level 1, or SL1 - Baird et al., 2009) microhabitats, each with a characteristic vegetation type, microtopographic relief, water-table and soil-moisture regimes, soil hydraulic properties, and soil biogeochemical regime (Ivanov, 1981;Alm et al., 1997;Belyea and Clymo, 2001;Rydin and Jeglum, 2006). These microhabitats commonly aggregate into larger spatial structures at horizontal scales of tens to 10 hundreds of m (SL2), often forming landscapes composed of strongly directional, nonrandom patterns that may be linear and contour-parallel, polygonal, or maze-like (e.g. Aber et al., 2002;Eppinga et al., 2008;Korpela et al., 2009). Baird et al. (2009) demonstrated that the frequency distribution of water-table depths from across a peatland landscape depends not just on the proportion of the landscape covered by different 15 microhabitats but also on the pattern. In consequence, pattern may play an important role in determining peatland-atmosphere fluxes of greenhouse carbon gases (carbon dioxide and methane), which vary with water-table depth (Bubier et al., 1993(Bubier et al., , 1995Roulet et al., 2007). Additionally, it seems likely that the development and maintenance of peatland patterning is governed by the same mechanisms that control peat accu-Introduction Observational (e.g. Foster and Fritz, 1987;Belyea and Clymo, 2001;Comas et al., 2005) and modelling (e.g. Nungesser, 2003;Swanson, 2007) studies have identified a variety of aspatial or one-dimensional (vertical only) feedback mechanisms that may help to explain directionless clumping of SL1 units into larger features. However, understanding the highly directional nature of patterning seen in many peatlands clearly 5 requires an explicit consideration of spatial interactions. This problem lends itself naturally to investigation using 2-or 3-dimensional simulation models, in which directional transfers of water, energy and nutrients, and their effects on pattern, can be explored directly. The idea of interplay between long-and short-range processes is a recurring theme in many apparently successful models of patterned landscapes, including peat-10 lands (Rietkerk et al., 2004a, b;Eppinga et al., 2009) and other landscape types such as marsh tussocks (van de Koppel and Crain, 2006) and dryland tiger bush (Lefever and Lejeune, 1997).
One current hypothesis on peatland pattern formation, the ponding mechanism, has been explored using cellular landscape models, but forms the subject of a curious 15 yet little-discussed disagreement between two groups of studies. The ponding mechanism consists of a feedback between water-table depth, peatland microhabitat succession, and peat hydraulic properties. Areas with deeper water tables are assumed to be more likely to support hummock vegetation, whereas areas with shallower water tables are assumed to be more likely to support hollow vegetation (cf. Rydin Eppinga et al. (2009) used their own cellular model to explore the individual and combined effects of the ponding mechanism and two other feedbacks. They found that contour-parallel striped patterning such as that shown in Fig. 1 developed under several combinations of these feedbacks, but not when the ponding mechanism alone was used; i.e. when using the ponding model in isolation they were unable to find a 5 combination of parameters that generated anything other than a homogeneous, unpatterned landscape, seemingly in disagreement with the earlier studies. Eppinga et al. (2009) employed in their models a scheme whereby every cell in the model landscape possesses continuous values of state variables such as hydraulic head, vascular plant biomass, and porewater nutrient concentration. This approach is distinct from the 10 binary hummock/hollow designation, and the discrete developmental timesteps during which microhabitat transitions occur, used in the SGCJ models. These and other differences in implementation between the SGCJ models and the model of Eppinga et al. (2009) mean that direct comparisons cannot strictly be made. However, the fact that Eppinga et al. (2009) were unable to generate striped patterning in their version of the 15 ponding model presents a curious disagreement between what would appear to be little more than alternative numerical implementations of the same essential conceptual model. Given that Eppinga et al. (2009) did not discuss this disagreement further, it is difficult to judge whether the ponding mechanism alone still represents a viable theory on peatland patterning.

Aim and objectives
We recreated a version of the ponding model to explore this disagreement. We investigated three characteristics of the SGCJ representation of the ponding mechanism that we suspected may have individually or collectively led to a Type-1 error in the SGCJ studies (i.e. causing a model to predict patterning despite the ponding mechanism be-Introduction  Table 1 for a glossary of algebraic terms), of any grid square in the original SGCJ models depends entirely on whether that square is currently a hummock or a hollow. Only two values of T are possible, one canonical value for hummocks, T hum , and one for hollows, T hol . A more realistic representation would have been to assign canonical 5 values of a more intrinsic property of peat such as saturated hydraulic conductivity, K hum and K hol [L T −1 ] to hummocks and hollows, respectively, and to calculate transmissivity as the product of K and the thickness of flow in that square, thereby allowing for continuous variation in T across the model landscape. We extended the hydrological submodel from the SGCJ models in this way in order to remove 10 any unrealistic constraints that the simplified hydrology may have placed upon overall model behaviour. We also wished to explore the effects of the absolute values of T or K (as appropriate) upon model behaviour.
ii. Impermeable deep peat: the original SGCJ models considered only a shallow layer of near-surface peat, and assumed that deeper peat was impermeable to 15 groundwater flow. The top few dm of peat are usually the most permeable (e.g. Fraser et al., 2001;Clymo, 2004) and are therefore prone to the most rapid subsurface flow, although deeper peat is rarely truly impermeable. Indeed, a number of studies have indicated that drainage through deep peat layers may play an important role in peatland development and the ability of these ecosystems to self-20 organise (e.g. Ingram, 1982;Belyea and Baird, 2006;Morris et al., 2011). The assumption of impermeable peat below the uppermost few dm may have prevented the SGCJ models from representing a potentially important hydrological interaction between surficial hydraulic structures and deeper peat layers. In our version of the ponding model we experimented with the effects of incorporating a 25 permeable lower layer.
iii. Dependence on hydrological steadiness: in the SGCJ implementation of the ponding mechanism, simulated water tables are, in theory, allowed to equilibrate Introduction with the current distribution of hummocks and hollows (and so the spatial arrangement of transmissivity) before the next microhabitat transition occurs (see Sect. 2 for full model description). The original SGCJ authors reported that the strength of patterning produced by their models increased strongly with the stringency of the criterion used to determine hydrological "steady state". The criterion considers the proportion of cells in which the rate of water-table change [L T −1 ] is less than a threshold rate. Steady-state conditions are deemed to occur when the proportion of cells below the threshold is greater than a proportion set by the model user. However, preliminary experiments (not reported here in full) with our own version of the SGCJ model indicated that the relationship between the hydrological 10 steady-state criterion and pattern strength may be more complex than previously reported. We chose to experiment with the effects of hydrological transience upon model patterning, either by manipulating the water-table steady-state criterion, or by driving the model with a real time series of rainfall data.

Overview
We began with a model that was as similar as possible to that employed by Couwenberg (2005), as far as the original model description allows. As well as this replica model (henceforth, Model 1) we created three additional models (Models 2, 3 and 4) with slightly altered routines in order to examine the effects of items i and ii, above (see 20 also Table 2). Each of the four models may be thought to consist of three submodels that simulate: shallow saturated groundwater movements and the spatial distribution of water- the spatial distribution of near-surface peat hydraulic properties. The new map of peat hydraulic properties is then used as an input to the hydrological model at the beginning of the next developmental step.

Hydrological submodel
The hydrological submodel uses a modification of the DigiBog model of peatland satu-10 rated hydrology . We took the original DigiBog Fortran 95 code and added routines to represent the ecological and hydrophysical submodels. In the current study we deactivated DigiBog's peat accumulation, decomposition and hydrophysical subroutines described by Morris et al. (2012). We refer the reader to  for a comprehensive description of DigiBog's governing equations and numerical im- 15 plementation, although a few points are pertinent here. DigiBog represents a peatland as a grid of vertical peat columns; in plan each column is equivalent to one square grid cell in the SGCJ models, although DigiBog also allows for vertical variation in peat hydraulic properties. Horizontal saturated groundwater flow occurs between adjacent columns (four-square neighbourhood) at a rate equal to the product of the harmonic 20 mean of inter-cell transmissivity (the product of depth-averaged hydraulic conductivity and the thickness of flow) and inter-cell hydraulic gradient, according to the Boussinesq equation (see also McWhorter and Sunada, 1977). The hydrological model is run for a predetermined length of simulated time in order to allow the water- to highly transient water tables that are still changing rapidly and are far from being in equilibrium with the hydrological inputs or boundary conditions; the opposite is true of long runtimes. In Models 1, 2 and 3 we varied ∆t e between 1 h (highly transient water tables) and 10 000 h (approximately 417 days; highly steady water tables) to examine the effect of the steady-state criterion upon model behaviour. We performed only a single simulation with Model 4, with ∆t e = 17 520 h (equal to two years of simulated time), and introduced hydrological transience by driving the model using a daily rainfall time series derived from observed data (see below).

Ecological submodel
In Models 1, 2 and 3 the probability, p [dimensionless], of any model cell being desig-10 nated as a hummock during a given developmental step is a linear function of watertable depth, Z final , at the end of the previous developmental step. When the water table is at the surface in any cell (i.e. when Z final = 0.00 m), p = 0 (i.e. that cell is necessarily designated as a hollow for the following developmental step). The value of p increases linearly with increasing water-table depth up to Z final = 0.05 m. When Z final is equal to 15 or greater than 0.05 m the cell is automatically designated a hummock (i.e. p = 1) (see Fig. 2). We refer to the upper 0.05 m of the soil profile as the transition zone. All cells that are not designated as hummocks are automatically designated as hollows. Swanson and Grigal (1988) provide a full description of microhabitat transitions as simulated in our Models 1, 2 and 3. In Model 4 microhabitat transitions are dealt with slightly dif-20 ferently. Rather than Z final , Model 4 assigns hummocks and hollows based on each grid square's time-averaged water-table depth, Z mean [L], during the second half of the previous 2 yr developmental step. The first 365 days of each developmental step are used to allow Model 4's simulated water tables to adjust to the newly updated distribution of peat hydraulic properties, but this adjustment period is not incorporated into Z mean . 25 This measure was intended to ensure that Z mean in Model 4 contains no artefact of water- water-table behaviour in an arguably more realistic manner than the short equilibration times used in Models 1 to 3.

Soil properties submodel
Hummocks are assumed to produce peat that is less permeable than that produced by hollows. DigiBog's implementation of the Boussinesq equation uses depth-averaged 5 saturated hydraulic conductivity, K [L T −1 ], and the thickness of flow to calculate saturated groundwater flux between adjacent columns. However, the original SGCJ models simply assigned a single value of transmissivity, T hum [L 2 T −1 ], to hummocks and another, higher value, T hol , to hollows. Model 1 uses the simple treatment of canonical transmissivities as per the original SGCJ models. We assumed default transmissiv- 10 ity values of T hol = 2.0 × 10 −4 m 2 s −1 and T hum = 1.0 × 10 −5 m 2 s −1 , thereby preserving the T hol to T hum ratio of 20 : 1 that led to strong patterning in the original SGCJ models. Models 2, 3 and 4 use a more sophisticated and realistic treatment allowed by DigiBog, whereby hummocks and hollows are assigned canonical values of hydraulic conductivity, K hum and K hol , respectively; T for each cell is recalculated during each iteration 15 of the hydrological submodel as the product of K and the thickness of flow, H (height of water table above the model's impermeable lower boundary) (cf. Freeze and Cherry, 1979). In this way T is able to vary in a continuous manner based on water-table position, which is more realistic than the simple binary-T treatment of the original SGCJ models. For Models 2, 3 and 4 we assumed default values of K hol = 1.0 × 10 −3 m s −1 20 and K hum = 5.0 × 10 −5 m s −1 , thereby giving a K hol to K hum ratio of 20 : 1. In models 1, 2 and 3 we also manipulated the default values of T or K by factors of between 0.05 and 20 in line with objective (i). In Model 1 we varied T hum between 5.0 × 10 −7 and 2.0 × 10 −4 m 2 s −1 ; and T hol between 1.0 × 10 −5 and 4.0 × 10 −3 m 2 s −1 , whilst always maintaining a T hol to T hum ratio of 20 : 1. Similarly, in Models 2 and 3 we varied K hum 25 between 2.5×10 −6 and 1.0×10 −3 m s −1 ; and K hol between 5×10 −5 and 2×10 −2 m s −1 , whilst maintaining a K hol to K hum ratio of 20 : 1. Introduction

Model spatial and temporal domains; boundary conditions
We implemented all simulations in a 70 (across-slope, x direction) × 200 (along-slope, y direction) grid of 1 × 1 m grid squares. In Models 1 and 2, the simulated peat aquifer overlaid a sloping impermeable base (replicating the assumption of an impermeable lower peat layer in the SGCJ models). Both the peat surface and the impermeable base 5 had a constant slope of 1 : 50; the permeable upper peat had a uniform thickness of 0.2 m. In Models 3 and 4 we assumed the thick, lower peat layer is also permeable with its own hydraulic conductivity, K deep , meaning that K for each cell is depth-averaged to account for the vertical transition in K between the upper and lower layers.  provide a full description of DigiBog's calculation of depth-averaged K and 10 inter-cell T . We used the groundwater mound equation (Ingram, 1982) to calculate the dimensions of a deep peat layer that is hemi-elliptical in cross section, has a uniform K of 1.25 × 10 −5 m s −1 , is 200 m from central axis to the margin, is underlain by a flat impermeable base, and receives a net water input (from an implied upper peat layer) of 155 mm yr −1 . The resulting deep peat layer was 3.97 m thick along its central axis, 15 curving elliptically down to the peatland's margin. Within DigiBog, this deep layer was overlain by a surficial, more permeable layer, which, like Models 1 and 2, was 0.2 m thick. We allowed all simulations to run for 100 developmental steps; the initial condition consisted of randomly-generated water tables, between 0.0 and 0.05 m below the peat 20 surface, in each cell. The initial microhabitat and soil-properties maps were based on this random initial water-table map in the usual manner described above.
We set the lateral (along-slope) boundaries to impermeable (von Neumann zero-flux condition) in order to simplify the representation of the simulated peatland within its local hydrogeological setting. In the simulations with an impermeable lower peat layer 25 (Models 1 and 2) both the upslope and downslope boundaries were set to constant water-table (Dirichlet) conditions of 0.025 m below the surface such that the overall hydraulic gradient along the model was equal to the topographic gradient and the gradient of the impermeable peat layer. As with the original SGCJ models, Models 1 and 2 received no rainfall; all inputs of water in these models came from the upper boundary condition, representing shallow groundwater influx and/or surface runon. In Models 3 and 4 we replaced the constant-head condition at the up-slope boundary with a no-flow (von-Neumann zero flow) condition to represent the drainage divide along the crest of 5 a raised bog in a manner similar to Morris et al. (2012) (see also Ingram, 1982), but retained the constant water- via rainfall to a simulated peatland, although it is pertinent to mention that U is assumed to be net of evapotranspiration; hence, its low value. We drove Model 4 using a daily time series of U. We took 365 days of daily rainfall data from close to Malham Tarn Moss, a raised bog in North Yorkshire, United Kingdom, for the calendar year 2011. The rain gauge recorded precipitation on 294 days of that year, with an annual 15 total of 1633 mm of precipitation. Maximum precipitation on any day was 78.6 mm on 10 August. We multiplied each daily rainfall total in the time series by approximately 0.245 so as to give an annual total of 400 mm of precipitation, equal to the constant U assumed in Model 3, whilst maintaining a plausible temporal variation (Fig. 3). During each 2 yr developmental step in Model 4 we cycled twice through the 365-day rainfall 20 series. By using the same year of rainfall data repeatedly in this way we were able to induce water-table transience without the potentially complicating effects of inter-annual variability in rainfall patterns.

Analysis
Previous studies (e.g. Andreasen et al., 2001;Eppinga et al., 2009) have demonstrated that the human eye is a powerful tool for assessing pattern strength in landscape models, and we used a visual appraisal of the two-dimensional microhabitat maps produced by the four models as our primary metric of simulated pattern strength. In addition we 42 Introduction  Swanson and Grigal (1988), as an objective and reproducible metric of pattern strength. The value of R increases with pattern strength: values less than 2 indicate unordered landscapes without discernible patterning; values of R greater than 4 represent highly ordered landscapes with clear, strong patterning.

5
We also calculated the proportion, S, of model cells occupied by hummocks during each developmental step. Finally, we calculated model turnover rate, Q, defined as the proportion of model cells that undergo a transition from hummock to hollow, or vice versa, in each developmental step.

Hydrological transience
In Models 1, 2 and 3, the strength of across-slope, striped patterning at SL2 initially increased strongly with increasing hydrological steady-state criterion, from an apparently entirely unordered, random mixture of hummocks and hollows at SL1 when ∆t e = 1 h, to peak "strengths" at around ∆t e = 100 h (for Models 1 and 2) or 50 h (Model 3). This 15 increase in patterning strength is evident from both a visual appraisal of the final microhabitat maps (Fig. 4), and an increase in the values of relative variance, R, with increasing ∆t e (Fig. 5a, b, c). For values of ∆t e greater than 50 h (Models 1 and 2) and 20 h (Model 3), an increase in the spatial scale of the simulated SL2 stripes is apparent with increasing ∆t e (i.e. the stripes became broader in the y (along-slope) direction).

20
Increasing ∆t e also led to SL2 stripes whose upslope and downslope edges were increasingly straight and sharply defined (Fig. 4). For values of ∆t e greater than 100 h (Models 1 and 2) or 50 h (Model 3), the simulated SL2 stripes had broadened to the point of appearing to be "over-developed", such that they no longer resembled realistic peatland striped patterning (cf. Fig. 1 As well as governing the spatial configuration of simulated patterning, the value of ∆t e also affected the temporal dynamics of Models 1, 2 and 3. For low values of ∆t e , turnover rates of model cells were high: when ∆t e = 1 h, approximately half of all cells would undergo a transition from designation as a hummock to a hollow (or vice-versa) during each developmental step, providing further indication of a random landscape.

5
Increasing ∆t e led to a reduction in turnover rates in Models 1 and 2 as SL2 structures began to stabilise, to a minimum of approximately 30 % of cells per developmental step when ∆t e = 100 h (Fig. 5d, e). Increasing ∆t e had little effect on turnover rates in Model 3 until ∆t e = 1000 or 10 000 h (Fig. 5f). The temporal trends in SL1 turnover rate and proportional hummock coverage also indicate what appears to be a limit cycle in the 10 behaviour of all models under highly stringent hydrological steady-state criteria. Particularly when ∆t e = 17 520 or 10 000 h, and to a lesser extent when ∆t e = 1000 h, the models exhibited a cyclical behaviour in time whereby the simulated landscape would alternate on an approximately regular cycle of between four and ten developmental steps between two contrasting states: a dry landscape with deep water tables, dom-15 inated by hummocks; and a wet landscape with water tables near the bog surface, dominated by hollows (Fig. 5). The contrast between the wet and dry states was particularly pronounced in Model 2, which cycled between approximately 15 % and 95 % hummock coverage when ∆t e = 10 000 h (Fig. 5h). The single simulation with Model 4 (∆t e = 17 520 h) behaved highly similarly to Models 1, 2 and 3 when ∆t e = 10 000 h, 20 predicting nearly homogeneous, unpatterned landscapes ( Fig. 6) that cycled rapidly between being dominated by hummocks and hollows (Fig. 5f, i).
In all simulations that generated striped SL2 patterning, those SL2 units migrated consistently downslope (not shown) in the same manner as that previously reported in the published accounts of the SGCJ models.

Simplified calculation of transmissivity
Any artefacts introduced to model behaviour by the simple, binary treatment of peat permeability employed in the original SGCJ models (see objective i, above) would have 44 Introduction been evident as differences in behaviour between Model 1 (simple, binary transmissivity scheme, based solely on hummock/hollow designation) and Model 2 (continuous treatment of transmissivity as product of peat hydraulic conductivity and thickness of flow). Models 1 and 2 behaved highly similarly to one another, in terms of their response to both changing values of ∆t e and changing absolute values of peat permeability. Both 5 models developed realistic looking, contour-parallel SL2 stripes over the entire model domain for intermediate value of ∆t e . In both models the striped patterning was weak and discontinuous for ∆t e = 15 and 20 h; the patterns were stronger and continuous when ∆t e was between 35 and 100 h; the SL2 stripes became unrealistically broad at ∆t e = 1000 h, and were mainly absent when ∆t e = 10 000 h (Fig. 4). Temporal patterns of summary metrics (SL1 turnover rate, hummock proportion, relative variance of hummocks per row) were qualitatively and quantitatively similar between models 1 and 2, and again responded similarly to changes in ∆t e (Fig. 5) and the ratio of T hum to T hol or K hum to K hol (not shown). 15 Any artefacts introduced to model behaviour by the simplifying assumptions made in the original SGCJ models about the hydrogeological setting of the simulated peatland (see objective ii, above) would have been evident as differences in behaviour between Model 2 (impermeable deep peat layer; constant-head upslope boundary condition; zero rainfall addition; constant slope of peatland surface) and Model 3 (permeable deep 20 peat layer; no-flow upslope boundary; constant net rainfall rate of 400 mm yr −1 ; hemielliptical aquifer shape). Model 2 simulations that generated patterning did so over the entire model domain, although the same was not true of Model 3. Striped patterning at SL2 developed in Model 3 only in the downslope portion of the model domain. Furthermore, the area of the model domain that exhibited patterning extended upslope approximately 140 m from the downslope boundary (Fig. 4). Clear, continuous striped patterning that extended all the way across the across-slope direction (x direction) of the model domain developed in Model 3 at much lower values of ∆t e than in Model 2. Continuous, closely-spaced SL2 stripes were evident in the downslope area of Model 3 when ∆t e = 20 h, although sharply defined, continuous SL2 stripes in Model 2 did not 5 develop for values of ∆t e less than 35 h. Additionally, Model 3 behaved differently from Model 2 in terms of its responses to changes in the absolute values of peat hydraulic conductivity (see below).

Manipulation of peat properties
Modest changes in the absolute values of T or K brought about large changes in 10 the nature and strength of simulated patterning in Models 1, 2 and 3, although the nature of these responses varied between models. In Models 1 and 2, the lowest values of T or K led to largely unpatterned landscapes, dominated by apparently random distributions of hummocks and hollows, with either a complete absence of patterning (×0.05 treatment) or very weak, discontinuous patterning (×0.1 treatment) (Fig. 7). In

Hydrological transience and ecological memory
The striking similarity between the patterns generated by our Models 1, 2 and 3 at intermediate values of ∆t e and the striped patterns reported by the previous SGCJ authors, as well as the absence of patterning in our models under more stringent steady-5 state hydrological conditions, are highly suggestive that the previous SGCJ authors had failed to attain a genuine steady state in their hydrological submodels. The dependence of realistic patterning in our models on intermediate equilibration times can be reconciled with the currently popular theory that landscape patterning generally arises from interplay between long-and short-range processes (e.g. Rietkerk and van 10 de Koppel, 2008; see also Turing, 1952). When equilibration times are very short each grid square's water-table behaviour is influenced only by its immediate neighbours, and long-range influences are unable to propagate across the model domain. Conversely, the longest (and arguably most realistic) equilibration times cause the simulated watertable map to be dominated by its long-range boundary conditions, which, under steady 15 state, eliminate much of the short-range effects of contrasting peat hydraulic properties in neighbouring SL1 and SL2 units. Without additional feedbacks such as those explored by Eppinga et al. (2009), it is only at intermediate values of ∆t e that the model strikes the balance of long-and short-range feedbacks seemingly required for patterning. 20 If the absolute values of ∆t e are taken literally then the model predicts that realistic patterning only occurs if micro-succession at SL1 operates on timescales of hours to days, rather than years. Such a prediction is clearly in error, meaning that it is perhaps initially tempting to discard the model in light of our work. However, we believe that this may be too rigid an interpretation of the hydrological steady-state criterion. It appears 25 that simulations with transient water tables have inadvertently captured some genuine aspect of peatland ecohydrology that leads to SL2 patterning. Moreover, Model 4, in which we introduced water- Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | through a daily rainfall series, did not generate patterning. This indicates clearly that it is the short equilibration times used in Models 1, 2 and 3, rather than simply nonsteady water tables (which are also present in Model 4) that caused patterning in our simulations.
For lower values of ∆t e , the water tables across the model's domain are not in full 5 equilibrium with the current distribution of hummocks and hollows. As such, those water-table maps reflect the distribution of hummocks and hollows during not only the current developmental step but also partly the previous developmental step. It may be that the formation of peatland patterning relies on some ecological memory effect (cf. Peterson, 2002), such as the influence of recently buried peat that no longer reflects 10 current positions of hummocks and hollows but which remains hydrologically important due to its shallow depth. Despite DigiBog allowing for 3-dimensional variation in peat properties, the ponding model is in essence a 2-dimensional model insofar as the new microhabitat map during each developmental step supersedes the previous map entirely; the modelled system retains no memory of previous SL1 units or their 15 associated soil properties. As such, our models neglect some of the structural complexity that peatlands exhibit in three spatial dimensions (e.g. Barber, 1981). Our simulations with non-conservative hydrological steady-state criteria (particularly ∆t e = 50 and 100 h) implicitly contain a type of ecological memory effect, because the influence of previous patterns of surface vegetation and shallow soil properties are expressed 20 in the model's hydrological behaviour; this memory effect decreases in strength with increasingly stringent hydrological steady-state criteria. Our results indicate that this memory effect is important to the formation of patterning in the ponding model and should be investigated directly in future studies, perhaps via the use of cohort-based peat accumulation models (e.g. Frolking et al., 2010;Morris et al., 2011Morris et al., , 2012. and ×20 treatments, the values of T or K were well within reported ranges; the fact that these parameterisations failed to produce patterning suggests that a full suite of relevant processes and feedbacks has not been represented. Particularly in Model 3, in which model water-table levels are maintained by the simulated addition of rainfall, even modest changes in the absolute values of K hol and K hum caused simulations to "run away" to either wet or dry end-member states. In simulations with higher values of hydraulic conductivity Model 3 drained rapidly to the downslope boundary, causing water tables to fall below the transition zone and giving rise to uniformly dry simulated landscapes dominated by hummocks. Conversely, the simulations with lower values of hydraulic conductivity caused Model 3 to drain so slowly that water tables rose to the surface of all columns, leading to uniformly wet landscapes dominated by hollows. The same effect was not evident in Models 1 and 2 because the constant-head condition at the upslope boundary maintained water-table levels within the transition zone. Nonetheless, manipulating T and K in Models 1 and 2, respectively, still resulted in landscapes devoid of realistic patterning, demonstrating high sensitivity to the absolute 15 values of those parameters. The high sensitivity to the absolute values of T and K may be taken to suggest that our models are lacking at least one important negative feedback, namely that between peat decomposition and changes in peat hydraulic properties. For example, saturated hydraulic conductivity is known to decrease strongly with increasing time-integrated 20 decomposition of peat (e.g. Boelter, 1969;1972;Grover and Baldock, 2013). Areas with low water tables are prone to more rapid decomposition and more rapid collapse of pore spaces. The resultant reduction in hydraulic conductivity in turn causes peat to retain water more readily and leads to reduced decay rates, stabilising system behaviour (see Belyea, 2009;. This concept is partially repre-Introduction  ). However, in the ponding model only two K values are possible, K hum and K hol , meaning that hydraulic conductivity cannot vary as a continuous function of decomposition; the ponding model's representation of this relationship may therefore be overly constrained. Modelling studies by Morris et al. (2011) and Swindles et al. (2012) have indicated that a continuous relationship between peat decomposi-5 tion and hydraulic conductivity may be highly important to the ability of peatlands to self-organise and to maintain homeostatic water-table behaviour. The representation of this negative feedback within a model of peatland patterning such as ours has the potential to stabilise model behaviour by allowing simulated water tables and peat permeability to self-organise, thereby reducing model sensitivity to small changes in soil 10 hydraulic parameter values. This would require the expansion of the ponding model so as to include routines that describe litter production and decomposition, and continuous changes in peat hydraulic conductivity. The inclusion of these processes in patterning models would also help to address the question of ecological memory effects raised above, and suggests the need to unify models of peatland surface patterning such as 15 those considered here, and cohort models of long-term peatland development (e.g. Frolking et al., 2010;Morris et al., 2011Morris et al., , 2012.

Improved model hydrology
Our alterations to the ponding model's hydrological basis compared to previously published versions had little impact on model behaviour, evidenced by the fact that Models 20 1, 2 and 3 all behaved in qualitatively similar manners, including their response to hydrological transience. The previous SGCJ authors reported that patterning is stronger and forms more readily as the slope angle of the model domain, and associated hydraulic gradients, are increased. The curved cross-sectional shape of the bog in Model 3 means that slope angle near the downslope boundary is greater than the uniform Model 3 then allowed patterning to spread upslope onto increasingly shallow slopes. The differences in behaviour between Models 2 and 3 can therefore be explained by the curved aquifer shape, reflecting the slope angle effect reported by previous authors. In particular, the nature of SL2 patterning in Model 3 (Fig. 4) is strikingly similar to that seen in the simulations with domed aquifers reported by Couwenberg and 5 Joosen (2005). We are therefore left to deduce that the permeable deep peat layer in Model 3 had little independent effect, and that the generation of patterning by the ponding mechanism is not dependent on either hydrological interaction with, or isolation from, deeper peat layers. The simple binary treatment of hummock and hollow transmissivity used in the earlier 10 SGCJ studies also appears to have produced no artefact in the behaviour of those models (or in our Model 1) compared to our more realistic Boussinesq treatment (in Models 2 and 3). Nonetheless, it is important to recognise that neither the simplified treatment of transmissivity nor the assumption of an impermeable deep peat layer, despite both being questionable assumptions in themselves, were responsible for the 15 model's reliance on the hydrological steady-state criterion and the absolute values of peat hydraulic properties, nor the prediction of downslope migration of SL2 stripes. As such we are confident that these behaviours (reliance on hydrological transience; high sensitivity to absolute values of transmissivity; downslope pattern migration) are genuine facets of the ponding model and are not artefacts of numerical implementation.

Pattern migration
The downslope movement of SL2 stripes in the ponding model is an emergent behaviour that should be treated as a testable hypothesis against which the ponding model could, in part, be tested (cf. Grimm et al., 2005). We are aware of very little direct evidence as to the long-term stationarity or otherwise of peatland patterning, 25 likely because the long timescales involved preclude direct observation. In one of the few studies into the matter Koutaniemi (1999)  evidence of consistent downslope migration of any of these features over a 21 yr period. Kettridge et al. (2012) used data from a combination of ground-penetrating radar and soil cores to investigate the below-ground physical properties of peat soil in a raised bog in Wales. They found subsurface layers that dipped consistently towards the bog's margins, which they interpreted as evidence of down-slope migration of SL1 or SL2 5 units during the bog's development. More observational work is clearly required to ascertain whether or not the findings of Kettridge et al. (2012) hold in the general case, but it appears that the ponding model's prediction of the downslope migration of SL2 stripes cannot be used to falsify the model at this stage. 10 Our numerical experiments suggest strongly that simulations reported in previous studies had not attained true hydrological steady-state conditions; under true steady-state conditions patterning does not occur. Realistic, striped patterns only form when microhabitat transitions occur in response to highly transient water-table patterns. The equilibration times required by the hydrological submodel to attain the level of tran- 15 sience necessary for patterning are of the order of hours to days: these timescales are largely meaningless in terms of vegetation dynamics. The models' reliance on the hydrological steady-state criterion indicates that ecological memory may play a role in pattern formation. Although ecological memory is not represented explicitly (or deliberately) in our models, the transient water tables appear to have inadvertently replicated 20 such an effect. We increased incrementally the sophistication and realism of the models' routines, chiefly by: representing a permeable deep peat layer; improving the models' calculation of transmissivity; and introducing a plausible temporal variation in rainfall. In doing so we removed a number of simplifying assumptions made in previous studies. How-Introduction sensitivity to hydrological steady state is a genuine facet of the models and is not an artefact of numerical implementation. The models appear to be unrealistically sensitive to the absolute values of the parameters used to represent peat permeability. We interpret this sensitivity as indicating that the models are missing a negative feedback between peat decomposition and hy-5 draulic conductivity, which previous studies have shown is important to the ability of peatlands to self-organise.

Summary and conclusions
Peatland structures and processes exhibit complexity in three spatial dimensions. A logical next step for patterning research would be to combine 2-dimensional (horizontal only) cellular landscape models such as those presented here with (mostly 1-10 dimensional; vertical only) peatland development models that provide a more detailed and realistic representation of peat formation, decomposition and dynamic changes in soil hydraulic properties. This would allow a direct process-based exploration of both the mechanism involved in ecological memory, and a continuous relationship between peat decomposition and hydraulic properties. Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | of discontinuous permafrost, Manitoba, Canada, Global Biogeochem. Cy., 9, 455-470, doi:10.1029/95GB02379, 1995 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Rietkerk, M. G., Dekker, S. C., de Ruiter, P. C., and van de Koppel, J.: Selforganized patchiness and catastrophic shifts in ecosystems, Science, 305, 1926Science, 305, -1929Science, 305, , doi:10.1126