Explicitly modelling microtopography in permafrost landscapes in a land-surface model (JULES vn5.4_microtopography)

Microtopography can be a key driver of heterogeneity in the ground thermal and hydrological regime of permafrost landscapes. In turn, this heterogeneity can influence plant communities, methane fluxes and the initiation of abrupt thaw processes. Here we have implemented a two-tile representation of microtopography in JULES (the Joint UK Land Environment Simulator), where tiles are representative of repeating patterns of elevation difference. We evaluate the model against available spatially 20 resolved observations at four sites, gauge the importance of explicitly representing microtopography for modelling methane emissions and quantify the relative importance of model processes and the model’s sensitivity its parameters. Tiles are coupled by lateral flows of water, heat and redistribution of snow. A surface water store is added to represent ponding

: Ice-wedge polygons at Samoylov (left, J. Boike) and palsas and mire at Iškoras (right, S.E. Chadburn). Both landforms form due to ground ice and show repeating patterns of elevation difference. Palsas have a larger elevations (0.5 -3 m) than polygon rims (up to 1m).

80
Microtopography can result in wetter areas in hollows which has a warming effect in summer through the increased thermal conductivity of saturated soil. Ponding can also have a warming effect by lowering the surface albedo. In winter wind-blown snow gathers in hollows and provides extra insulation. The combined effect is that lower areas tend to be warmer and wetter, which in turn results in greater emissions of methane, and permafrost thaw in those areas.

85
Microtopographic features are highly prevalent in the permafrost region, as they can be generated as a result of small-scale processes through the formation or removal of ground ice. While our methods are generally applicable, this study applies our model to two such features, ice-wedge polygons and palsas ( Figure 1). Ice wedge polygons are created where wedge ice forms in cracks created by the thermal contraction of the ground, pushing up the adjacent ground (Lachenbruch, 1962). Polygons are widespread; wedge ice is thought to be ~20% or more of the uppermost permafrost volume (Liljedahl et al., 2016), and are 90 found in the continuous permafrost zone. Palsas are formed where a localised cold patch of ground causes frost heave, elevating mounds of peat (Seppälä, 1986). These in turn experience shallower snow cover and further heaving. When these mounds are joined to cover a larger area, they are known as peat plateaus. Palsas / peat plateaus are found in the discontinuous and sporadic permafrost zones, and are also widespread, though not to the extent of polygons Kremenetski et al., 2003;Fewster et al., 2020). While it is hard to get an estimate for the total extent of palsas and peat plateaus, permafrost peatlands 95 are estimated to cover 1.7 million km 2 . The choice of these two landforms affords the opportunity to test our model in climatically distinct regions and using different configurations.
Microtopography can preserve permafrost, but, conversely, can cause localised areas of raised ground temperature or 'hotspots' where abrupt-thaw events can be initiated. These hotspots can experience slumping due to ground ice thaw, extending and 100 exacerbating the hotspot and possibly forming or extending a waterbody -a thermokarst lake (Langer et al., 2016). These abrupt thaw events expose new carbon to decomposition and could have a large contribution to the permafrost carbon feedback, potentially doubling the radiative forcing from carbon fluxes as result of gradual thaw (Walter Anthony et al., 2018). Such abrupt thaw features are common, with thermokarst lakes covering an estimated 11% of Siberian Yedoma and much of Northern Canada (Walter et al., 2006). However, even rarer types of abrupt thaw events are hypothesised to have a large 105 impact. Thaw-related hillslope erosional features are projected to affect a tiny 3% of abrupt thaw terrain, but be responsible for a third of abrupt thaw emissions . Explicitly modelling microtopography could therefore be important as a prerequisite to the inclusion of abrupt thaw features and their effect on the permafrost carbon feedback in land surface models.

110
Modelling efforts have been hampered by a lack of long-term high-resolution observations. Nevertheless, at a handful of field sites the local variability has been well quantified and at these sites 'tiling' approaches to modelling micro-and mesotopography have been tested as a route to representing sub-grid processes in ESMs (Langer et al., 2016;Aas et al., 2016 of permafrost landscapes, for which they have shown that including micro-and even meso-topography is essential. When 115 applied on the global scale, a more top-down rates-based approach, similar to that taken by Schneider von Deimling et al., 2015) may well be necessary due to the complexity of the feedbacks involved. However, this results in a loss of flexibility and predictive power, which could result in the omission of potentially important feedbacks. We therefore believe that an explicit tiling approach should be further pursued and developed. In this paper, we therefore follow the approach of Aas et al. (2019) and Nitzbon et al. (2019) and construct a scheme within JULES where two interacting tiles, offset by an 120 elevation difference, are taken as representative of a repeating microtopographic unit. The relative areas within this repeating unit are then representative of the entire grid-cell ( Figure 2). The two tiles interact via exchange of both soil moisture and heat, by surface run-on and by redistribution of snow. A surface water store is added to the low tile to represent ponding. The twotile approach was chosen as the simplest configuration that still includes an explicit representation of microtopography, as our focus is testing this first step-up from the previous implicit microtopography of JULES. For the same reason, dynamically 125 changing areas and elevations, necessary for modelling abrupt thaw processes, as well as a thorough treatment of waterbodies, though important, are not included in the scope of this investigation. Rather, this work aims to strengthen the foundations for applying tiling approaches to microtopography in permafrost landscapes in three main ways. Firstly, by evaluating the efficacy of a two-tile model against microtopographically-resolved observations of snow depth, soil temperature and moisture at four sites. Secondly, by evaluating the importance of explicitly modelling microtopography with regards to its effect on the 130 modelled methane emissions. Thirdly, by quantifying the effect of individual model processes, the model's sensitivity to individual model parameters and the resulting uncertainty in the modelled summer soil temperature as a result of uncertainty in those parameters. This sensitivity testing will enable future modellers to assess which model processes should be focused on, which could be simplified, and which parameters most need to be constrained. 135 Figure 2: The implicit microtopography of standard JULES (left) vs our explicit approach (right). In standard JULES (left), the gridbox average values of soil temperature and carbon are used to calculate a methane flux. This is then multiplied by a wetland fraction, which is calculated by TOPMODEL using the modelled water table depth and the gridcell topographic index. Our explicit approach (right) uses two soil tiles offset by an elevation difference to represent repeating patterns of microtopography. Model spatial parameters are shown in grey and are defined for individual sites in Table 2. 140 2 Methods

Sites
Four sites were used for model validation: two Siberian polygon sites, Samoylov and Kytalyk, and two Scandinavian palsa sites, Stordalen and Iškoras. This choice facilitates comparison, as sites for each type of landform have similar climate, while polygon sites are colder and drier than palsa sites. Table 1 shows a summary of climatic variables for each site and further site 145 information can be found in Table A 1 (Appendix A). Data used from each site can be found in Figure A 1 and Table A 2 (Appendix A). Samoylov, Kytalyk and Stordalen have already been set-up and validated in JULES (Chadburn et al., 2017) and all of these sites have microtopographically resolved data. https://doi.org/10.5194/gmd-2021-285 Preprint. Discussion started: 26 October 2021 c Author(s) 2021. CC BY 4.0 License. Samoylov and Kytalyk are sites with ice-wedge polygons, where observations of soil moisture and temperature used in this paper are for low-centred polygons. Samoylov island is situated in the Lena river delta . Measurements are from the Holocene terrace on the Eastern part of the island. The soil is characterized by alluvial deposits (65%), with the 155 organic layer underlain by sand / silt with some peat layers. A shallow gradient in the terrace leads to a gradient in hydrological conditions, and consequently variation in the degradation of polygons (Nitzbon et al., 2020a). Areas of high and low-centered polygons can be seen, along with a large areal fraction of water surfaces (25%, Muster et al., 2012). Kytalyk is in a region of Yedoma, with Pleistocene ice deposits occurring nearby in 20 to 30 m tall hills (van der . The measurement site itself is in a depression which was originally a Holocene thermokarst lake. 2 to 3 m below this is the current river plain, 160 where active thermokarst features are present alongside the river. In addition to ice-wedge polygons, Kytalyk has areas of small palsa-like hills . The elevation differences caused by frost-heave have led to general difference between dryer elevated areas and lower flooded areas, however ponding is only seen in small areas (up to 10%) where ice wedges are actively developing. Precipitation is smaller than that of the palsa sites, but evaporation is limited to the short growing season, meaning the soil can remain wet. 165 Iškoras and Stordalen are sporadic permafrost sites with palsas surrounded by a lower-elevation, non-permafrost mire. Both have been experiencing lateral degradation of palsas. Stordalen has experienced warming of 2. 5°C between 19135°C between and 20065°C between (Callaghan et al., 2010, and has lost up to 10% of the palsa between 1970 and 2000 (Malmer et al., 2005). Stordalen has areas of palsa, ombotrophic bog which receives runoff from the surrounding palsa, and fen. Drainage out of the palsa follows the 170 frost table via narrow 1 to 2 m channels (Olefeldt and Roulet, 2012). Unlike the palsa and bog, the fen has no permafrost. In addition to receiving drainage from the bog, the fen receives a large amount of water from a lake just outside the peatland complex, and therefore the 1.7km 2 catchment surrounding it. The Iškoras peat plateau is small, covering an area of ~4 ha, and is surrounded by and interspersed with fen and ponds. Similar drainage channels are observed in palsas, along with typical signs of palsa degradation such as cracking and blocks of peat detaching from palsa edges . 175

JULES
JULES (jules.jchmr.org) is a community land surface model (in  and Clark et al. (2011)). It is the land surface scheme used in the UK Earth System Model (UKESM) (Sellar et al., 2019) and can also be used as an offline model, both globally and at the site level as is done here. For this study, version 5.4 (http://jules-lsm.github.io/vn5.4/, last access 2 nd June 2020) was modified. In JULES a gridcell has a single soil column, on which surface tiles representing different plant functional 180 types and / or non-vegetation tiles such as bare soil, water and land ice can be used. Each surface tile models its own multilayer snowpack and canopy water store, and fluxes of heat, water, carbon and nitrogen are aggregated before exchange with soil layers. For this study plant competition is modelled by the TRIFFID dynamic vegetation scheme. While the standard configuration of JULES usually has four soil layers with a depth of 3m, here we use a configuration with 20 soil layers with a total soil depth of 7.8 m as suggested by Chadburn et al. (2015). JULES has been used to model permafrost at the site level, 185 as well as its global extent and the permafrost carbon feedback, including wetland and methane emissions (Chadburn et al., 2015;Burke et al., 2017b;Comyn-Platt et al., 2018;Burke et al., 2017a;Gedney et al., 2019).
The version of JULES used here includes two tiles representing either the palsa and mire or the polygon rim and center. This is done by repeating the same driving data twice within the model. The following developments were included within JULES 190 to improve the representation of microtopography in the northern high latitudes. We have added snow redistribution and lateral water and heat fluxes between the tiles and a surface water store to represent ponding on the lower tile. In addition, we describe the JULESspecific modifications for this model, which are a more conceptually consistent parameterisation of the water table depth diagnostic and a new method of limiting layer numerical over-or under-saturation. Model parameters will be defined as they are encountered but can also be found in Table 2. 195

Snow redistribution
To model the effects of snow blown by wind on the snow depth of each tile we use the scheme of Aas et al. (2019) andNitzbon et al. (2019) where snowfall, ( −2 −1 ) is instantaneously partitioned between the two tiles according to the current snow depth on the tiles. We assume that while the current snow depth, ( ) on either tile is less than the snow catch, catch = 200 0.05 , the snow is 'caught' by vegetation, so cannot be redistributed. In this case, each tile experiences the same snowfall.
Snow above this depth is not caught, so if the level of snow on the adjacent tile is lower, snowfall will be instantly redistributed to it. A buffer, buffer = 0.03 is applied, such that redistribution does not happen if the elevations differ by less than this amount. This means that once the snow is approximately level, snow is once more able to fall on both tiles. This recreates the observed behaviour of snow filling hollows first (Wainwright et al., 2017). Finally, as the areas of the tiles may be different, 205 and snowfall in the model is per area, the redistributed snow is subject to an area scaling factor, 1 2 , where 1 and 2 ( 2 ) are the areas of each tile. The full scheme is therefore: High to low tile redistribution: Low to high tile redistribution: where Δ ( ) is the elevation difference between tiles, and the subscripts 'high' and 'low' specify the tile with which the 210 variable is associated.

Lateral water fluxes
Lateral flows of water were introduced into JULES using an approach mirroring the existing calculation of vertical fluxes, where fluxes are calculated based on the difference in matric potential between soil layers due to their level of saturation.
Existing functions for calculating the layer matric potentials and hydraulic conductivities are used, and the fluxes interfaced 215 into the existing code for the water balance for each layer. Unlike for the vertical fluxes, no implicit correction is used, as fluxes are assumed to be relatively small. Here, fluxes are calculated for horizontal connections, where connections occur between a layer and any other layers horizontally adjacent to it, taking into account the area of overlap of each connection ( Figure 3). However, a 'sloped' flow scheme where layers are sequentially connected to their corresponding layer in the neighbouring tile was also considered. Further discussion on the strengths and weaknesses of the two schemes can be found 220 in the supplementary material. In JULES, the vertical water flux between layers per land surface area, is calculated using Darcy's law and is of the form: where ℎ ( −2 −1 ) is the hydraulic conductivity, Δ ( ) is the distance between the centres of the layers, and Δ ( ) the difference in matric potential between layers. When flows are horizontal, the matric potential gradient instead uses the 225 horizontal distance, Δ ( ), and as the gravitational potential, ( ), remains constant, Δ Δ = 0 and the 1 can be dropped.
As with the snow scheme, we must go from the intensive form used by JULES of fluxes per land surface area, to the extensive form of −1 through the interface, and back again. The calculated flux is therefore multiplied by area of the interface, Δ Δ , to go from −2 −1 to −1 , where Δ ( ) is the vertical overlap (orange in Figure 3), and Δ ( ) is the contact length. This divided by the area of the tile that is being considered, 1 or 2 , so that the change in water content 230 of a layer is still expressed per unit area of that tile. The areal scaling factor for tile = 1 or 2, , is therefore: where ( 2 ) denotes the area of tile . The lateral flow for a connection per surface area of tile ℎ , ( −2 −1 ) is then calculated as:

235
In JULES, and ℎ can be calculated using either the Brooks and Corey (BC) or the Van Genuchten (VG) relations (see . Here we use the BC relation for ℎ . The VG relation is used for , as it avoids a mismatched when both layers are saturated but is otherwise almost identical to the BC relation for the soil properties of the modelled sites. As in the case for the vertical fluxes, the conductivity is determined using the average liquid water content of the connected layers. This calls into question the use of only two boxes, as the sloped area of polygonal tundra has been observed to have a 240 lower saturation than either the polygon rim or centre  which could lead to the flow being restricted. Not modelling the slope is however similar to previous approaches Nitzbon et al., 2019) and is the simplest configuration that still models microtopography.
The horizontal flow scheme raises the question of what to do with the upper layers of the raised tile, which may have no 245 horizontally adjacent layer. If the water table in the elevated tile is above the surface of the lower tile, and above the level of the surface of any ponded water (if present), water will be able to laterally egress the soil (Figure 3 C). Similarly, if a pond is present on the low tile and the surface of the pond is above the level of the water table in the high tile, water can flow from the pond laterally into the soil (Figure 3 D). For these flows it is therefore necessary to determine the level of the local water table.
For the high tile, for a particular connection, the local pressure head in the high tile ℎ ( ) is the height of the water table  250 above the midpoint of the vertical overlap (orange in Figure 3 is no matric suction from air. For water ingress, Δ is the height of the pond above the midpoint of the connection, plus the matric suction of the layer the water is entering. While the local pressure head must be calculated in order for horizontal flows between the soil and the air or the pond and the soil to be modelled, flows are always from a saturated or partially saturated region to an unsaturated one. The model therefore does not implement saturated lateral flows, which are conceptually important for flows between polygons (Wales 265 et al., 2020). This means that advective flows of heat may not be properly represented. However, the model will still be able to act to balance the water table and moisture potentials between tiles, though the rate with which equilibrium is reached may be different. Any discrepancy in rate is however expected to be less than the usual timescales over which the water table changes and will therefore not be a problem. The addition of saturated-saturated flows would require solving for the pressure potentials of multiple saturated layers of varying soil properties, which may or may not contain frozen water. This is non-270 trivial, as the pressure potentials themselves depend on the (now 2D) water fluxes, leading to the requirement for a more complex iterative scheme. On balance, we consider this scheme to be a reasonable simplification for our purposes.

Lateral heat fluxes
Lateral fluxes of heat by conduction and advection use the same connections as the water fluxes. However, of these, only the soil to soil connections are used, as the pond thermodynamics are not explicitly modelled (it is assumed to be at the surface 275 temperature). The equations for lateral heat flows have the same form as those for the vertical fluxes between layers, as before with the addition of the appropriate areal scaling factors, namely: where ( −2 ) and ( −2 ) are the conductive and advective heat flows for a lateral connection per surface area of 280 the tile respectively. As before, refers to the tile that is being considered.
is the lateral water flow from the previous section, Δ ( ) is the temperature difference between laterally adjacent layers and ( −1 −1 ) is the specific heat capacity of water. The average thermal conductivity of the two soil layers ( −2 −1 ), is calculated using the usual JULES function , using the mean of the frozen and unfrozen water contents and of the dry soil thermal conductivity for the two layers being considered. 285 In reality, the heat flux is determined by the continuous derivative / . The temperature profile could vary horizontally in a different way to the moisture profile, and indeed would be expected to be asymmetrical across a frozen-unfrozen interface.
In the future, the need for Δ could be eliminated by dynamically determining this profile and including lateral freezing/thawing at the interface. For now, in the model freezing and thawing are effectively distributed evenly throughout the 290 layer undergoing phase change and we use the same Δ as in the water flux calculation. https://doi.org/10.5194/gmd-2021-285 Preprint. Discussion started: 26 October 2021 c Author(s) 2021. CC BY 4.0 License.

Ponding and run-on
In the studied sites, the lower area tends to remain almost saturated year-round ( Figure 6 & Figure 8). However, in standard JULES the surface layers rarely become saturated, as after rainfall events or snowmelt water quickly infiltrates into a space provided by evaporation, or, as happens after snowmelt, water is in excess of the maximum infiltration rate and runs off. At 295 both types of site, ponding is common Klaminder et al., 2008), both in the centres of polygons and around the edges of palsas. This enables water to gradually infiltrate rather than running off and provides a buffer that enables the soil surface to remain saturated. To model this, a surface water store was added to the low tile. At each timestep, a fraction, , of runoff from the high tile joins the throughfall of the low tile as run-on, along with the amount of water in the surface water store. Water in excess of the maximum infiltration is then stored in the surface water store. When a pond is present, bare soil 300 evaporation is switched off, and evaporation from the surface water store is equal to the potential evaporation rate calculated by JULES . Currently, no thermal effects of the pond are included, and the pond cannot freeze, though in future aspects of FLake (Rooney and Jones, 2010) could be used to introduce these processes, in a similar manner to (Langer et al., 2016). The ability to pond improves the surface wetness, and if the surface is able to become saturated year-round, the pond can build up to become a persistent waterbody ( Figure 6). When not limited, even in the wettest configuration, pond 305 depths did not rise to more than a couple of metres. However, as this may not physically possible due to the size of the hollow, ponded water in excess of a fraction pd of the elevation of the high tile is removed as runoff.

Diagnosing the water table position in a partially saturated layer
The implemented lateral flow scheme requires the position of the local water table for a layer to be determined, here referred 310 to as the 'layer head', or . Here, we present a new method of determining the water table within a partially saturated layer in a way that is consistent with the assumptions of the cell-centred numerical method for solving the Richards equation used by JULES. The model procedure for finding the layer head is then described. As mentioned in the previous section, JULES calculates the matric potentials of layers using the fractional saturation of the layer, . Previously, the position of the water table in JULES has been calculated by assuming that refers to the saturation at the midpoint of the layer, . This aligns 315 with the usual calculation of vertical Darcy fluxes using the potential difference between the midpoint of each layer. Setting the potential of the water table to be 0, the soil moisture above the water table is assumed to have the equilibrium profile such that the increase in gravitational potential is balanced by a decreasing matric potential such that the total potential is constant.
Using the Brookes-Corey (BC) relation at the layer midpoint, where ( ) is the height above the water table, ( ) is the BC matric potential at saturation and is the Clapp and Hornberger (1978) soil exponent. The result of this calculation is shown by the blue curve in Figure 4. A similar calculation can be done using the Van Genuchten (VG) relation, for which we use an alternative formulation: with similar results to BC, giving the orange curve. This correctly places the water table at the midpoint when saturated, rather 325 than being offset by .

Figure 4:
The fractional saturation of the layer, , vs the distance of the water table below the top of the layer − . Previously, the position of the water table has been calculated by assuming that refers to the fractional saturation at the midpoint of the layer, resulting in the blue and orange curves, and the water table not exceeding the midpoint of the layer 330 when the layer is saturated. However, in JULES, when considering how much water fits in a layer, refers to the average saturation of the layer. Using this interpretation results in the green curve, where the water table is at the top of the layer when the layer is saturated. Making the simplification that the green curve is approximately linear while the water table is within the layer, the position of the water table within a partially saturated layer can easily be found from the average saturation of the layer. 335 However, in JULES, is also a prognostic variable tracking the average saturation of a layer and is updated according to how much water flows into or out of a layer. This contradicts the calculation of the water table depth and of the direct calculation of the matric potential at the midpoint from , as a saturated layer ( = 1) would have its water table at the top, and conversely the saturation at the midpoint could be 1 even when is not. We therefore need a way of calculating the water table from the value of the average saturation of the layer. We start by doing the reverse. If the water table distance below the top of the layer, 340 − (m), is known, the average fractional saturation of the layer, , can be found by dividing the integration of the equilibrium VG profile for over the layer by the layer thickness, ( ), noting that below the water table is saturated.
This gives the green curve in Figure 4. The water table is now at the top of the layer when the layer is saturated. The deviation from either the standard VG or BC calculation is small when the water table is below the bottom boundary of the layer, indicating that ≈ is a good approximation for most cases away from the saturated/unsaturated boundary. Accounting for the difference is therefore only of any importance in wetland environments. To invert the function to calculate the water 350 when < : The final method is simple to implement and compares well to the full numerical integration ('VG integrated' in Figure 4), 360 with fully accurate end-point values.
Now that the position of the water table within a partially saturated layer can be found, the model procedure for finding the local water table for a layer is as follows: if a layer is saturated, then each layer above the original layer is checked sequentially until a layer that is partially saturated is found. The position of the water table within this partially saturated layer is then 365 determined, and the total height from the original layer to this is the layer head. If the original layer is partially saturated, the position of the water table within the original layer is determined only if the layer beneath is saturated. If a layer has frozen water present, then the layer head is zero for that layer. If, while sequentially checking the layers above the original layer, a layer which contains frozen water is found, the layer head is the height from the original layer to the bottom of that layer. If the saturation continues to the surface, then the layer head is the height from the original layer to the pond surface (if present). 370

Limiting over-or under-saturation
The standard methods JULES uses to avoid supersaturation or undersaturation as a result of the water flux calculation numerics can result in the unintended consequence of water being passed out of the soil column. This is particularly a problem for freezing saturated soils. Here, we present a method which avoids this problem, and which also integrates with the scheme for simulating lateral fluxes of water. 375 In some cases, the explicitly and implicitly calculated Darcy fluxes could lead to layers which have a fractional saturation outside the range of 0 to 1. To avoid this, JULES can be set to either limit the water flux into the top of the layer (l_soil_sat_down = false) or to limit the flux out of the bottom of the layer (l_soil_sat_down = true) so that the saturation does not go outside this range . The soil water extraction and sub-surface runoff are not limited, as they cannot 380 cause oversaturation, and should be absent in the case of under-saturation. However, both approaches can cause problems, as in both cases the choice of which flux to limit does not account for the direction of water flow. This can cause water to either be able to pass downwards through permafrost, or be ejected upwards out of the soil when water is drawn upwards towards a frozen surface layer when refreezing (Figure S 1 and Figure S 3, Supplementary). In addition, an approach is needed that can also limit the lateral fluxes when necessary. 385 Here, we solve this problem by only the limiting the incoming fluxes when a layer will potentially be oversaturated (soilsatupdown). This is done by multiplying them by an appropriate factor such that the layer becomes exactly saturated (Figure S 2, Supplementary). A multiplicative factor was used so that all fluxes scale equally. For a potentially under-saturated layer, any outgoing fluxes (discounting qbase and extraction) are limited in a similar manner. This is necessarily an iterative process as 390 for any layer either flux could be limited, potentially creating the need to adjust a flux in an adjacent layer. The iteration was found to be most efficient sweeping upwards from the lowest layer, and repeats until either there is no appreciable change, or the iteration limit (set at 1 greater than the number of layers) is reached. Generally, only a few of sweeps are required, and the iteration limit is never reached. As lateral flows can be limited, which changes the flows of the paired tile, the limiting code is run a second time for each box after it has been run for all boxes once. 395 https://doi.org/10.5194/gmd-2021-285 Preprint. Discussion started: 26 October 2021 c Author(s) 2021. CC BY 4.0 License.
For the most part, soilsat-updown results in a soil moisture profile very similar to that of l_soil_sat_down = false, and there are only a very few cases where flux out of the top of the soil is avoided ( Figure S 1, Supplementary). This does however suggest that l_soil_sat_down = false is in general the more physically realistic scheme. This being said, the few instances where the top of the soil is saturated, and layers of soil exhibit an upward suction into saturated layers are of importance to 400 freezing wetland environments and become more frequent with the inclusion of ponding.  point out that, considered globally, l_soil_sat_down = false generates too much runoff in permafrost environments due to correctly allowing a perched water table above a frozen layer. However, the perching of the water table above saturated permafrost is a key behaviour responsible for permafrost wetlands (Avis et al., 2011), so setting l_soil_sat_down = true is undesirable. Our addition of ponding will help reduce runoff, while the effect of macropores (Bechtold et al., 2019) should also be considered, 405 so as to get the right infiltration in the right places for the right reasons. As a side note, an alternative, more physically realistic solution to water being drawn into already saturated frozen layers would be the inclusion of frost heave.

Site parameters
The model is parameterised using observed characteristic dimensions of landscape features present at sites. These parameters 410 are: the elevation difference between tiles Δ ( ), the areas of each tile 1 and 2 ( 2 ), the contact length between tiles Δ ( ) and the horizontal distance between tiles Δ ( ) ( Figure 5). For the polygon sites, the raised and lower tiles represent the polygon rim and centre. For the palsa sites, the tiles represent the palsa and mire. The ratios of these parameters determine the model's geometry, so the model could be readily applied to other forms of microtopography also. We have approximated ice-wedge polygons to have a circular centre surrounded by a ring, so 1 = 1 ( Δ 2 ) 2 , while palsa and mire are represented by 415 squares of equal area with one side in contact, so 1 = 2 = Δ 2 . Δ is chosen to be representative of the length over which the transition between the thermal and hydrological regimes of the two tiles takes place, as it is this that affects the magnitude of the fluxes. Δ is therefore less than the distance between the centres of each tile.

420
Polygon spatial parameters for Samoylov were found by measuring polygons in the Digital Elevation Map by Julia . Averages were taken over the major and minor axes of each polygon, and these used to calculate an average area for the rim and centre. For Kytalyk, the polygon Lhc11 in (Teltewskoi et al., 2019) was used. In both cases, these are low centred polygons. For the palsa-mire configuration, the geometry of two bordering squares was chosen as a simple 425 configuration representing the edge of a peat plateau, this being an area of interest due to being the site of palsa degradation.
In reality, the geometries of peat plateaus are complex, leading to extended perimeters. Indeed, both the perimeter to area ratio and the rate of areal degradation have been observed to increase as palsas become smaller (Borge et al., 2016). The outcome of varying Δ will therefore be of interest for these sites. The chosen geometry for peat plateaus can however be a valid representation if a number of assumptions are true. That the border is straight is a good approximation if the border length is 430 small. That the tiles only exchange fluxes on one side relies on gradients of moisture and temperature being small on the other https://doi.org/10.5194/gmd-2021-285 Preprint. Discussion started: 26 October 2021 c Author(s) 2021. CC BY 4.0 License. three sides. This ought to be true of the two sides perpendicular to the border due to symmetry, and the border being small.
Gradients of soil moisture and temperature will also be small at the side furthest from the border if the distance between this side and the border is great enough. Squares of side length 5 m were chosen for both sites as a compromise between these assumptions. The model's sensitivity to this choice was tested later. The elevation used for Iškoras was taken to be the same 435 as the elevation of the palsa where the snow depths were measured (0.68 m) so that results could be validated. This falls towards the lower end of the range for the site observed by Martin et al. (2019), namely 0.1 to 3 m, but is nonetheless an established palsa. Stordalen is a site experiencing degradation, for which Karlgård (2008) gives the characteristic height of 0.5 m, Olefeldt and Roulet (2012) give a range of 0.5 to 2 m, Klaminder et al. (2008) a range of 1 to 3 m. Due to the lack of specific microtopographic observations for this site, an elevation difference of 2 m was chosen to provide a contrast to the 440 model setup for Iškoras. The resulting estimates of site-specific model parameters are given in Table 2.

Configurations of JULES simulations
For model code and configuration files, see the code and data availability statement. Each site was spun up in standard JULES (without any microtopographic additions or tiling) for 10,000 years using repeating driving data for the years 1901-1910.
Admittedly, this length of spinup turned out to be unnecessary as methane fluxes were ultimately calculated offline using observed carbon profiles. This spun up state with the additional modification of a saturated soil column was used to initialise 450 the main run (1901 -2016) for both standard and tiled JULES. As the same spun up state was used to investigate multiple configurations which may have different hydrological end states, the soil column was set to be initially saturated for both the spinup and the transient run. This was to avoid the possible development of dry layers during spinup persisting due to being frozen, which could then limit the possible end states (Figure S 5 and Figure S 6, Supplementary). The thermal and hydrological states of the tiles of the transient runs were stable by 1911, and our analysis is limited to the 21 st century. Forcing data were 455 prepared as described in (Chadburn et al., 2015), where reanalysis data for the grid cell containing each site were adjusted using site observations where available. Water and global change Forcing Data (WFD) were used for the period 1901-1979 (Weedon et al., 2010(Weedon et al., , 2011 and WATCH Forcing Data Era-Interim (WFDEI) used for the period 1979-2017 (Weedon et al., 2014(Weedon et al., , 2018. Parameters for the main tiled runs are given in Table 2 & Table 3. For the sensitivity study, three different approaches were used to test the size of the effect of different model processes and the model's sensitivity to its parameters: • Additive/subtractive process switching: a series of runs where aspects of the model were switched on individually, using standard JULES as the base configuration, or switched off individually using the full tiled runs as the base 465 configuration.
• Latin hypercube sampling: a series of 100 runs where model parameters were varied independently using a Latin hypercube sampling, each parameter having an even distribution between the maximum and minimum values given in the columns S.S. min and S.S. max in Table 2. These runs enable the possible range for an output variable to be gauged, with the caveat that the low ratio of samples to parameters mean that configurations where a number of 470 parameters are at their extreme are unlikely to be sampled.
• Individually varying parameters: A series of runs where we varied each parameter individually using same maximum and minimum values as for the hypercube sampling, while holding the other parameters at their standard tiled configuration values (as in Table 2). These runs enabled the effects of individual parameters to be disentangled. Using these outputs and estimates of site-specific parameter uncertainties enabled the uncertainty in an output variable due 475 to the uncertainty in a particular parameter to be found.

'Offline' methane analysis
The methane fluxes for each site were estimated 'offline' in R using the equations from the methane model in JULES applied to soil carbon and temperature data, similar to part of the analysis in S. E. Chadburn et al. (2020). The difference here was that we used a hybrid approach where observed soil carbon profiles were used in combination with modelled soil temperatures (as 480 opposed to using only observations). We took this approach as it removes error introduced via poor simulation of soil carbon, since the soil carbon was not developed or evaluated in this study. This approach enables the impact of substrate vs temperature dynamics to be disentangled. For this work we look at methane emissions as calculated by the previous non-microbial version of JULES, as well as by the recently added microbial scheme. The microbial scheme adds a dissolved organic carbon pool, a methanogen pool, and the associated dynamics of methanogenic growth and dormancy. This has been shown to better represent 485 the observed seasonal dynamics (S. E. Chadburn et al. 2020). Work to improve the soil carbon model in JULES is ongoing.

qbase configuration (modelling baseflow)
JULES uses TOPMODEL to calculate the baseflow (qbase) or 'subsurface runoff' for each soil layer, based on a topographic index and the position of the water table. This acts as the gridcell-scale lateral drainage. Here, we also qbase set to zero for layers which are unsaturated and/or frozen. However, in this explicit representation of microtopography, the lateral flows 490 between tiles directly model at least part of this flow. For this reason, we include the parameters 1 and 2 , which scale the drainage calculated by TOPMODEL applied to the high tile and low tile respectively. We first considered the configuration where qbase is set to zero for the simulated palsa and for the polygon centre and is on for the other tiles. This configuration is referred to as 'part qbase', or 'qbase on' as this is the main configuration where we allow qbase in the model (Table 3). The motivation for this is that there cannot be any drainage from the centres of polygons without going through the rim, and 495 similarly drainage from palsas would usually pass through the surrounding mire. There are two main problems with this approach. Firstly, we are applying a model based on gridcell-scale assumptions to a subgrid model. This may result in qbase being too great or too small. For instance, rather than being calculated using the gridbox average saturation, it is now calculated using the saturation of either the wetter or the drier part of the landscape. This approach could also result in a reduced qbase as it is only being applied to part of the area. Secondly, the lateral drainage calculated by TOPMODEL may not be suitable for 500 some instances of wetland, or where a large proportion of the gridcell is wetland. This is partly due to TOPMODEL's https://doi.org/10.5194/gmd-2021-285 Preprint. Discussion started: 26 October 2021 c Author(s) 2021. CC BY 4.0 License. limitations in flat landscapes (Koster et al., 2000), but is also due to wetlands being formed in areas where lateral drainage is impeded, or even which are fed by groundwater. In JULES, qbase only ever removes water, so the latter is not possible. The repeating nature of landscape units provides motivation for some degree of impedance, as two symmetrical bordering units can have no flow between them. As we expect that qbase may be further impeded, we also consider the configuration, referred 505 to as 'no qbase', where qbase is switched off for both tiles. This is not a general solution and is a departure from standard configurations of JULES. In reality, the lateral landscape-scale drainage will depend on the connectivity of the polygon troughs or mire network (Liljedahl et al., 2016;Connon et al., 2014), in addition to the presence of external reservoirs and inter-gridcell flows. These things are beyond the scope of this study and will require further consideration. In addition to these two main configurations for qbase, the response of the model to varying fd1 and fd2 between 0 and 1 was tested in the sensitivity study 510 (as described in section 2.5.2 Configurations of JULES simulations), and a configuration with qbase on for both tiles ('twin qbase') was included in the model process switching tests.

Results and discussion
An overview of the model output can be seen in Figure 6, which displays climatologies for Samoylov and Iškoras, showing soil moisture and temperature with depth, along with snow depth and modelled pond height, comparing observations with 520 standard and tiled JULES. These results will be expanded in the following sections, which examine the effect of the tiled model on snow, before moving on to soil moisture and temperature, and considering which model processes have most effect on these variables. The resulting effect on methane is considered, in addition to possible applications of the model. Finally, the sensitivity of the model to its parameters is investigated, and sources of uncertainty identified. Where climatologies are shown, the years shown are for those the model is averaged over, while the years averaged over for observations will be for all years 525 available (see Figure A 1  Iškoras, comparing observations with standard and tiled (no qbase) JULES. The white line above the plots of soil temperature shows simulated snow depth, the cyan line above the plots of simulated soil moisture shows simulated pond depth. For Samoylov, the tiled model has little effect on thaw depths. For Iskoras, thaw depths are much greater, and while the low tile does not get as warm as observed, the tiled simulation clearly differentiates between the permafrost palsa and non-permafrost mire. Iskoras is a wetter site, and the tiled simulation is able to develop a persistent pond. The greater elevation difference at Iskoras means that the saturated low tile has less effect on the upper layers of the high tile. The tiled simulation does not develop a pond for Samoylov, however, while not shown on this figure, a seasonal pond does develop for the other polygon site, Kytalyk. Soil moisture for the low tile is improved for Iskoras vs standard JULES, however the simulation for Samoylov remains too dry.

Snow depths
The splitting in snow depths between the tiles as a result of the snow redistribution scheme is shown in Figure 7. When comparing the average climatologies for Iškoras, where microtopographically resolved timeseries of snow depth are available, both the snow depths for the top and the base of the palsa are over 20 cm closer to observations for March than the standard 530 JULES simulation. Snow depths for the polygon centre for Samoylov are similarly closer. The goodness of the model's operation mostly depends on the forcing applied, how well the typical maximum snow height on the elevated areas matches the snow catch, and the tile areas. This is because most of the time the modelled low tile snow depth for these sites does not reach the elevation difference, the point at which snow in the model would build up on both tiles at the same rate. Consequently, the snow depth on the high tile is usually limited by the snow catch parameter ( ℎ ), while the low tile snowfall receives a 535 fixed adjustment of snow for the time that the snow depth is above the snow catch. An alternate 'sloped' snow scheme was trialled whereby the fraction of snowfall redistributed varied linearly between all the snowfall being redistributed if the https://doi.org/10.5194/gmd-2021-285 Preprint. Discussion started: 26 October 2021 c Author(s) 2021. CC BY 4.0 License. difference in snow elevations equalled the elevation difference between tiles, to no snow redistributed if the snow elevations were equal. The motivation for this scheme was that it might have been more suitable for polygon sites where the slope of the rim gradually becomes more covered in snow as the depth of snow in the centre increases. However, there was little difference 540 between average snow depths for each scheme, so the 'sloped' scheme was found to be unnecessary. Some of this discrepancy may be due to the forcing data used. For instance, in the forcing data for Samoylov there is almost no snowfall prior to the start of October, whereas observations indicate snowfall before this, resulting in around 5cm too little depth of snow for October to January. The scheme also has some obvious omissions. While snow is assumed to be redistributed by wind, wind speed and direction is not input into the snow scheme. Spatial variability due to wind direction and drifting 555 against slopes is therefore not considered. Much of the snow can also be blown a greater distance and accumulates in larger depressions such as lake basins and river gullies, therefore in reality the total balance of snow needs to be considered on a much larger scale. However, a more complete accounting for the spatial variability due to effect of wind speed and direction, topography and vegetation, as in SnowModel (Liston and Elder, 2006), moves towards needing a detailed gridded representation of topography to run, and goes against our aim for a simple representation of subgrid processes. Currently, the 560 snow catch is unaffected by the type and height of vegetation present in the model. Linking the snow catch and vegetation height would be a key next step and would enable snow-vegetation feedbacks to be investigated. Lastly, the properties of the falling snow have no effect on the modelled redistribution, and no consideration is given to wind erosion of the snowpack at a later time after the snow has settled. Zweigel et al., (2021) included the properties of falling snow and the effects of internal snow processes and snow microstructure on snow mobility for a high-Arctic site on Svalbard. They found that a thinner layer 565 of lower density snow built up on their high 'ridge' tile than their lower 'snowbed' tile, which was then able to be subject to wind erosion at a later time. However, while including these omissions would affect the properties, timing and spatial variability of the modelled snow redistribution, the snow depths on each tile have still been successfully differentiated with https://doi.org/10.5194/gmd-2021-285 Preprint. Discussion started: 26 October 2021 c Author(s) 2021. CC BY 4.0 License. the current model. Collection of additional spatially resolved time series data would enable further validation, investigation of the effect of missing processes and for model parameters to be better constrained. 570 Figure 8 shows in black the observed average saturation profiles for July, where sites are grouped according to landscape type due to a lack of detailed observations of hydrology at Kytalyk and Stordalen. For all studied sites, the lower area is characterised by year-round saturated or almost-saturated conditions, with a water table generally within 10 cm of the surface. Polygon rims are still fairly wet, with the observed water table (shown by box plots) within 20 cm of the surface, while observations from 575 palsas have an unfrozen liquid fractional saturation of around 0.6 at 10-20 cm depth. Model outputs in Figure 8 will be discussed in the following sections.

The importance of impeding baseflow (qbase)
Following the discussion on the treatment of baseflow in section 2.5.4 qbase configuration (modelling baseflow), model results in Figure 8 are shown both with the TOPMODEL calculated qbase switched on in the top two panels and switched off in the 580 lower panels. This switch needs to be considered before examining the other aspects of hydrology because of the potential size of its effect, the uncertainty in its implementation, and its departure from the standard hydrology of JULES. Switching off qbase has a large effect on the palsa sites, but a small effect on the polygonal sites, where the thaw depth is shallow ( Figure   8). We remind the reader that we have made the modification that qbase is always only permitted from saturated, unfrozen layers. For the palsa sites and in standard JULES, impeding qbase is enough to significantly raise the water table, saturating 585 the soil at 20 cm depth for Iškoras. However, the soil only becomes fully saturated with the tiling scheme (for the low tile, with qbase switched off). Conversely, if qbase is not impeded for these sites, the tiling scheme has negligible effect on the saturation profile. Impedance of qbase is therefore a necessary condition for modelling palsa mire sites correctly in JULES.
This agrees with the simulations of Martin et al. (2019), who when varying the external water flux in their peat plateau model found a "drainage effect" when transitioning from an external influx of +1.5 mm/day (+550 mm/year) to an outflux of -2 590 mm/day (-730 mm/year). This was enough to change the soil from saturated to well drained, decreasing the 1 m deep mean annual temperature by up to 2°. With qbase on, qbase in the tiled scheme becomes over -200 mm/year on average for Stordalen and over -500 mm/year on average for Iškoras ( Figure B 3, Appendix B), which is at the drained end of the transition found by Martin et al. We originally hypothesised that the calculated qbase could be too high for the low tile for palsa mires, due to being calculated using the saturation of the wetter tile (section 2.5.4 qbase configuration (modelling baseflow)). While 595 qbase is significantly larger from the low tile than from standard JULES (~500 vs ~150 mm/year respectively for Iškoras), qbase is still high enough to create a drained landscape even in standard JULES. From the sensitivity study, the individually varying parameters runs showed that by varying 2 (the fractional scaling applied to the calculated qbase for the low tile) the point at which saturation is reached can be seen (Figure C 1A, Appendix C). In fact, very little drainage is required for the mires to become unsaturated, and as such there are very few parameter configurations from the Latin hypercube sampling runs 600 that achieve saturation for these sites. The sensitivity of the model to the calculated qbase points to the need for more accurate determination of baseflow from permafrost wetlands. For the rest of this paper, the tiled configuration with qbase switched off ('Full tiling (no qbase)') is regarded as the best performing model configuration and is the one presented in subsequent   . Indicative observed water table depths for palsas are from measurements of 8 hollow points, not including ponded areas, and 5 palsa points in September 2007 . Water table depths for July from hollows in Iškoras were available but are not shown (C.T. Christiansen, H. Lee, unpublished data), as the range of values is contained within the range shown for Stordalen. The no 610 qbase tiled runs represent the saturation profiles at palsa sites well, with water table depths within the range of observations. The same is true for the low tile for Kytalyk. Switching off qbase appears necessary for the low tile for palsa sites but has little effect on polygon sites. Lateral flows out of the palsa site high tile have a similar effect to that of the TOPMODEL calculated qbase. Observed soil moisture for the palsa at Iškoras at 0.4 m decreases vs that at 0.2 m due to still being partially frozen in July. Note that simulations show the average saturation for the layer whose midpoint is at the depth shown.

The effect of microtopography on soil moisture and fluxes
Compared with standard JULES with qbase off, the tiling scheme causes a splitting in the level of relative saturation (fraction of saturation) (Figure 8). For polygons, the raised tile remains broadly unchanged compared with standard JULES, whereas the lower tile gets wetter. For the palsas, the lower tile becomes slightly wetter and can now be fully saturated. The raised tile becomes much drier due to being able to laterally drain onto the lower tile, and the larger elevation difference of the palsa sites 620 mean that the soil layers are less affected by the saturation of the low tile or even the level of the pond. Lateral drainage from the raised palsa tile has a very similar effect to that of qbase when switched on from the low tile or in the standard JULES simulation. In both cases, the lower tile is improved while the high tile remains, or becomes, too dry. Altogether, the palsa sites seem fairly well represented, though as we have seen, most of the credit for this is due to the disparity in drainage. While  Compared to standard JULES, the high tile drying appears to be driven by lower snow mass causing decreased snowmelt water influx, as well as the addition of lateral egress out of the soil. Fluxes are balanced by decreased runoff due to the decreased snowmelt. The low tile wetting appears to be driven by runoff decreasing almost to zero and the addition of run-on, an increase 635 in snowmelt water influx due to increased snow mass, and the lateral egress of water from the high tile to the pond. The change in these fluxes is balanced by large amounts of pond evaporation. For Stordalen, there is also a small amount of lateral flow through the soil from the low tile to the high tile, though this is negligible for other sites. Changes in runoff and snowmelt water influx, due to spatial variability in snow depth as well as run-on and ponding, appear therefore to be the main drivers of model soil moisture heterogeneity in this model, with lateral egress from the soil playing a smaller role at equilibrium. 640 Unsaturated soil-to soil flows are negligible, due to a combination of horizontally adjacent soil layers being saturated and/or frozen. While lateral unsaturated fluxes may have greater relative importance if the landscape drains or becomes drier, the hydrological difference between tiles would be less pronounced in this case anyway. This suggests that simpler, water-table based models of lateral flow may be adequate in wetland environments.

Distinguishing the impacts of model processes on soil moisture 645
In order to directly attribute the effect of different model processes on soil relative saturation, a series of runs were performed with individual processes switched off in turn with tiled no qbase as the base configuration, the 'subtractive process switching' https://doi.org/10.5194/gmd-2021-285 Preprint. Discussion started: 26 October 2021 c Author(s) 2021. CC BY 4.0 License.
runs. The resulting effect on the relative saturation is shown in Figure 10. We find that although individually the effect of each process is small (aside from the effects of qbase, previously discussed, and lateral flow from the high tile in palsa mire sites, dark blue in Figure 10) most processes have an effect in some conditions, and also that the importance of different fluxes varies 650 with the degree of saturation. For example, if the mire is well drained or already fully saturated, ponding will have little effect.
For polygons, switching off snow redistribution has one of the largest effects, and reduces the low tile liquid water content as a fraction of saturation, henceforward referred to as fractional saturation, by ~− 0.05 to −0.1. Polygon hydrology is also affected by the thermal state of the ground. For Samoylov, switching off lateral flows of water and heat (conductive and 655 advective) has the largest effect, reducing the splitting in level of fractional saturation from 0.2 to 0.06. While this is the result of switching off two flows rather than one, the size of the effect is still surprising when considering that turning off lateral water flow alone has little effect and that the effect is larger than the reduction in splitting for switching off lateral heat alone.
Indeed, this is the case even though the average lateral fluxes of water in the standard tiled run are negligible, and the corresponding temperature splitting is very small (Figure B 3, Figure 11). Similarly, for Kytalyk, switching off lateral heat 660 fluxes has the next largest effect after snow redistribution. This is perhaps unsurprising considering that this effect singlehandedly reduces the temperature splitting from 2° to almost nothing. Switching off runoff and ponding have little effect.
However, when the opposite method is tried and aspects of the model are switched on individually using std JULES as the base configuration (the 'additive process switching' runs), ponding and runoff can have a larger effect ( Figure B 1

The effect of microtopography on soil temperatures
The two types of microtopography have different effects on soil temperatures. Figure 12 shows the observed and modelled 675 soil temperatures for July at 0.19 m. The model summer temperature splitting is smaller than observed for palsas (3.5 vs 5°C), with the mean remaining lower and approximately unchanged compared with standard JULES, while polygons display small (0 to 0.5°C) temperature splitting, in agreement with observations. Thaw depths are likewise similar for high and low tiles at polygon sites, while palsa sites show a clear difference between palsa and mire (simulating permafrost and non-permafrost conditions respectively). In the model, Iškoras shows signs of being on the edge of palsa degradation, with a talik forming for 680 three consecutive winters after a warm summer in 2012 and greater winter snowfall for those years. Figure 11 displays climatologies for soil temperatures at 0.19 m depth, showing that the model is able to reproduce the overall behaviour for both types of landscape over the course of a year. Year round, the observed differences in temperatures between the polygon rim and centre are comparatively small (< 2. 5°C at 20cm depth). However, polygon rims are observed to get 685 colder in winter, be quicker to warm up in summer, can be warmer in summer, and cool down more quickly in winter. This general behaviour is reproduced in the tiled model, though the model lacks a substantial zero-curtain in the autumn when refreezing. This is probably due to the modelled polygon sites being too dry. Only the low tile of Kytalyk approaches having a substantial zero-curtain behaviour, it also being the polygon tile nearest to the observed saturation. Microtopography tends to https://doi.org/10.5194/gmd-2021-285 Preprint. Discussion started: 26 October 2021 c Author(s) 2021. CC BY 4.0 License. split temperatures so that they fall either side of that of standard JULES. If the temperature modelled by standard JULES lies 690 outside of the observed temperature splitting and is greater or smaller, the temperatures modelled by the tiled run also tend to be biased in the same direction and greater or smaller respectively. While polygon sites are close to observations in July, both standard JULES and the tiled configuration are generally colder than observed in summer, in addition to warming up and cooling down earlier. In addition to the aforementioned differences in soil moisture, part of the discrepancy in timing may be due to the shorter period of forcing data snowfall. That Samoylov is too cold in winter may be in part due to having too little 695 snow. Palsa mires show a more pronounced splitting in temperature for both the summer and the winter than polygons, and the model does well at capturing the difference in winter temperatures for Iškoras, and the mire zero curtain. However, the temperatures for both standard JULES and the tiled runs are too small in summer. The tiled model reproduces the seasonal temperature splitting behaviour, though lacks a substantial zero curtain for the polygon sites due to being too dry, and summer temperatures are cooler than observed for all sites. Soil temperatures for the tiled runs tend to lie either side of those for standard JULES. Figure 11 shows the effect on soil temperatures of the additive process switching runs, where model processes are switched on one at a time starting from the base configuration of standard JULES. Turning on the snow redistribution only results in almost the full winter temperature splitting for palsa sites, and too great a temperature splitting for polygon sites. The effect of snow redistribution on summer temperatures is smaller, though it is still responsible for the majority of the temperature splitting. While their effect is smaller than that of snow, drainage, ponding and lateral flows of water increase the summer 710 temperatures of the low tile, with the order of their respective magnitudes of effect approximately following that of their effects on soil moisture (Figure 12, Figure B 1 Appendix B). In particular we can see the difference between the best estimate and best estimate (no qbase) runs for the palsa mire sites, showing that the full temperature splitting cannot be achieved without a saturated low tile, though it remains to be seen if the high tile would also get significantly warmer if it too were wetter. https://doi.org/10.5194/gmd-2021-285 Preprint. Discussion started: 26 October 2021 c Author(s) 2021. CC BY 4.0 License. While its effect is not visible in Figure 11, lateral heat conduction has the effect of reducing the temperature splitting by -0.6 715 to -1.1° for palsa sites, and -1.1 to -2.0° for polygon sites (Figure B 2, Appendix B). This eliminates almost the entire temperature splitting for Kytalyk. There are almost no advective heat fluxes (Figure B 4, Appendix), however this may be due to the lack of saturated soil water fluxes in our model. The surface sensible and latent heat fluxes for the high and low tiles remain approximately the same as for standard Jules (Figure B 4, Appendix B). The modelled lateral heat conduction per unit length of border for Iškoras is 2.9 Wm -1 . Over a year, the energy transferred per metre of border could melt 0.28 m 3 of ice. 720

Distinguishing the impacts of model processes on soil temperatures and quantifying heat fluxes 705
Interestingly, this is of the same order of magnitude as the 0.13 m 3 yr -1 m -1 of volumetric loss of permafrost peat plateau edge found by Martin et al., (2021). Of course, not all the lateral heat flux would go into melting ice, so this rate of thaw would be smaller. Nevertheless, it suggests the possibility that the modelled heat flux between just two tiles could be used to calculate the rate of lateral thaw at the interface, and hence the rate of areal change from one tile to the other. by the boxes. The model goes some way towards reproducing the observed temperature difference between tiles, though the modelled temperatures and temperature splitting are smaller than observed. Out of the model processes shown, snow redistribution and switching off qbase have the largest effect for palsa sites. For Kytalyk the effect of ponding and lateral flows can also be seen. The effect of lateral flows of heat cannot be seen in this figure, as it only has an effect where some other factor is already causing a temperature difference. 735 Figure 13 shows methane fluxes calculated using the observed soil carbon profiles combined with soil temperatures from the 'standard JULES' and 'tiled (no qbase)' runs. In this figure, standard JULES has qbase switched on while tiled runs have qbase off. This is to represent the difference in methane fluxes between how JULES is usually configured and a fully microtopographic approach. Based on the effect of the tiled model on temperatures seen in the previous section, we might 740 expect to see a larger difference in methane production compared to standard JULES for the mire, and little difference with https://doi.org/10.5194/gmd-2021-285 Preprint. Discussion started: 26 October 2021 c Author(s) 2021. CC BY 4.0 License. respect to the polygon centres. Indeed, the total methane production for the polygons in JULES remains virtually unchanged for Samoylov and only 8 or 9% higher for Kytalyk. The palsa sites show an increase in methane fluxes vs standard JULES of 9 or 10% using the non-microbial scheme (NM), however when using the microbial scheme (M) for Iškoras this increases to 49%. This is likely a factor of whether or not the microbes can survive/grow in all soil layersfor example they are less likely 745 to grow a population in permafrost soil layers. These results show that explicitly modelling microtopography can increase modelled methane fluxes vs the standard JULES in some cases. As previously mentioned, here we have used the observed, rather than modelled soil carbon profiles to calculate the methane flux. In reality, the carbon substrate might also be greater in the wet parts of the landscape compared to the landscape average (Wagner et al., 2003;Eckhardt et al., 2019), so the effect on simulated methane that we show here may be amplified when soil carbon dynamics are modelled. 750

760
When comparing the simulated fluxes with flux tower observations (dotted lines in Figure 13), fluxes show reasonable agreement for palsa sites, though slightly lower than observed for Samoylov, but are more significantly lower than observed for Stordalen. This is what could be expected due to average modelled low tile summer temperatures for July and August being < 1° lower than observed for polygons, but > 4° lower for palsa sites (Figure 11). Previous modelling using a standard version of JULES with the modification of a saturated soil column has been able to better reproduce the observed flux tower 765 methane emissions for Stordalen (Chadburn et al., 2020, Figure 3). This underlines the primary importance of the core model setup in accurately modelling methane emissions, and that the difference made by modelling microtopography is only a refinement to this. In this case, the difference may be that the very upper layers of the lower tile in the model do not remain saturated year-round ( Figure 8) unlike the fully saturated column in Chadburn et al., (2020). This drier surface layer of soil provides an insulating effect, meaning that the soil does not heat up as much in summer. For the full tiling (no qbase) runs, 770 Stordelen's pond is only seasonal, developing after snowmelt and disappearing before midsummer. In contrast, Iškoras has a persistent pond and a fully saturated lower tile, resulting in higher summer temperatures and methane fluxes more comparable to those observed for Stordalen.
While eddy covariance methane observations (EC obs.) have the benefit of providing a detailed timeseries and a landscape 775 aggregate, they do not resolve small scale variations. For Figure 13, in order to compare the EC obs. with the modelled fluxes for the lower tile, for Kytalyk and Stordalen the EC obs. were selected for when the tower footprint was over the wetter area of the landscape, while for Samoylov the fluxes were divided by the observed wetland fraction (Table A 2, Appendix).
Chamber measurements are able to resolve these spatial variations, however it is harder to get a landscape-scale average and measurements often lack an extended timeseries. The average simulated fluxes are within the ranges of chamber flux 780 measurements for Samoylov (2.7 ± 1.3 2 ℎ −1 , ) 1 ) and for Kytalyk (3 ± 2 2 ℎ −1 , (Huissteden et al., 2005) 2 ) though these are not shown on Figure 13 due to the lack of a multi-year average. Methane chamber measurements from Bäckstrand et al. (2010) are shown for Stordalen. These measurements are the average from 2002-2007 for the 'green season' (days 119-288). Chamber measurements have been grouped as in (Johansson et al., 2006;Bäckstrand et al., 2010), where the semiwet area (three chambers) corresponds to Sphagnum spp. and Carex spp. and the wet area (two chambers) to 785 Eriophorum angsifolium. Three chambers were used for the palsa microsite. The wet area is completely saturated and permafrost free. The semiwet area has a water table fluctuating within 20 cm of the surface, and while this area not raised like the palsa, it is only partially thawed. A large difference in methane emissions can be seen between the wet and the semiwet areas. Our modelled lower tile fluxes are hard to directly compare with these areas. The lower tile sits in between the two in terms of moisture and arguably in terms of temperature (based on thaw depths), and so the methane fluxes could be 790 approximately correct for an area of the mire which is between wet and semiwet. However, as the saturated areas contribute a much larger share of the total methane emissions (see e.g. orange vs blue horizontal lines on Figure 13D), simulating methane emissions from unsaturated parts of the landscape will have little impact on the upscaled fluxes.
Another problem for upscaling is that the method of calculating methane emissions used here is now decoupled from changes 795 in the water table. In reality, the thickness of the oxic layer determines the amount of substrate unavailable to methanogens 1 Three low-centred polygon centres of three chambers each, uncertainty shows standard error in spatial variability. Seasonal average (10 th July to 18 th September 2006). This study confirmed that the chamber measurements matched up to those measured by the eddy covariance tower using the wetland fraction used here. 2 Two low-centred polygon centres (27 th -30 th July, 2004). Large variability can be seen for fluxes from the surrounding area.
https://doi.org/10.5194/gmd-2021-285 Preprint. Discussion started: 26 October 2021 c Author(s) 2021. CC BY 4.0 License. and affects how much of the generated methane is oxidised. This is represented in standard JULES by varying the wetland fraction depending on the water table level, as described in the introduction. Here, the low tile area is fixed, and so our calculation of the fluxes will increasingly be an overestimate as the water level drops below the surface. Future versions of the model will therefore require changes in the water table and the methane fluxes to be linked. As Chadburn et al. (2020) observe, 800 doing so would enable the microbial dynamics to be better modelled, as methanogens would be unable to survive in layers where the saturation is varying rapidly. Adding water table dynamics would also add to the feedback on the methane fluxes described previously of an insulating dry layer limiting summer soil temperatures, making the low tile emissions more sensitive to changes in the water table. Again, this is a contributing factor to the observed spatial variability in methane emissions and must be considered when upscaling. 805

Sensitivity study
The sensitivity study investigated the response of the model to varying model parameters using the methods described in section 2.5.2 Configurations of JULES simulations. Using the estimated possible ranges of parameter values, we determined the range of possible modelled soil temperatures for each site and quantified how the uncertainty in model parameters affects the uncertainty in model output. These tests are important as they confirm that the model will not give unreasonable results if 810 the parameters are set differently, while also indicating which parameters are in most need of constraining. Figure 11 and Figure 12 show the variability in the temperature splitting (TS, the difference between high and low tile temperatures) for July at 0.19 m as a result of varying parameters using Latin Hypercube Sampling (LHS). It can be seen that this variability is comparable to the size of the observed TS itself. For the polygonal sites, the choice in parameters can result 815 in the TS being larger by up to ~3°, as the usual splitting is small (0.5° for Samoylov). A similar increase from the usual modelled splitting is also possible for the palsa sites, however as the modelled TS is already smaller than observed, the resulting splitting can only be up to ~1.5° larger than observed. The parameter ranges used for the LHS are generous, and so we can be confident that including microtopography is 'safe' and will not result in unreasonable modelled temperatures. However, that the variability in the TS is similar to the size of the modelled TS means that care must be taken when choosing parameter 820 values. Table 4 gives approximate estimates of parameter uncertainties by site, which were then used alongside the results from the individually varying parameters runs (as in Figure C 1C, Appendix) to find the resulting effects of the parameter uncertainties on the size of the July TS at 0.19 m depth. The parameter uncertainties are estimates of our uncertainty in the parameters used, 825 and not a quantification of site level variability based on any large-scale survey. For Kytalyk, our choice of parameters appears to have resulted in a local minimum for the TS, such that for four of the parameters, variation in either direction results in an increase of magnitude of the TS. As a result, either tile could become the warmer tile, so microtopography could either increase or decrease the fluxes. The largest resulting change in the TS as a result of the uncertainty in an individual variable is 1.5° for the palsa sites and 0.7° for the polygon sites. This can be a large fraction of the overall TS and is larger than the effect of 830 many of the individual model processes, reinforcing the need to carefully chose model parameters. This is particularly the case for 1 (the high tile area, which affects the depth of snow on the low tile), 1 (fraction of drainage from the low tile, which has been previously discussed), and Δ (the horizontal distance, which is hard to measure but effects the magnitude of lateral fluxes). At the other end of the scale, uncertainties in the values of (maximum pond depth), 1 (drainage from high tile), and (run-on fraction) have little individual effect (≤ 0.1°) on the TS. There can be a large variability in palsa heights at 835 an individual site, for instance, Olefeldt and Roulet (2012) record a range of 0.5 to 2 m for Stordalen. It had therefore been a concern that a single pair of tiles may not be representative if this also corresponds to a large variability in temperature.
However, Figure C 1C, (Appendix) shows that there is little change in the modelled temperatures once the palsa elevation is https://doi.org/10.5194/gmd-2021-285 Preprint. Discussion started: 26 October 2021 c Author(s) 2021. CC BY 4.0 License. above 0.5m, meaning that having only two tiles can be a valid representation of multiple palsas whose heights are above this threshold. While we have taken care in choosing model parameters based on observations, it is recommended parameter 840 validation be implemented for future users, so that unreasonable or even unphysical combinations of parameters are not allowed. The sensitivity study also provides further motivation for a dynamic representation of the parameter Δ , as this currently responsible for a 0.6° uncertainty in the TS.  The results of the individually varying parameters runs also give insight into the effect of varying parameters on winter soil 850 temperatures, and particularly the model response to parameters affecting the snow scheme. While the main behaviours follow directly from the model equations (1 & (2), it is worth noting that soil temperatures are most affected by changes in snow depth when snow depths are shallower. For instance, soil temperatures are highly sensitive to the snow catch parameter ( ) for the high tile, where the sensitivity tests show that a difference of 5 cm can result in as much as 2° difference at 0.19 m ( Figure C 1.B, Appendix). However, the same is not true of the response of the low tile temperatures to snow catch as the snow is deeper, 855 so changes have less effect. Similarly, changes to 1 (the area of box 1) are a key control on snow depths on the low tile and can affect the snow depth by a much greater amount but effect the temperature less. For instance, changing the 1 by 50 m 2 for Stordalen results in changing the low tile snow depth by around 0.5 m and a temperature increase of around 1°.

Summary and outlook
In this study we have constructed and tested an explicit representation of microtopography within the JULES land surface 860 model. We have followed the two-tile approach suggested by Aas et al., (2019) and Nitzbon et al., (2019), though in order to implement such a model in JULES we have taken a different approach to implementing lateral fluxes, runoff and ponding. Our purpose has been to test the efficacy of an explicit representation of microtopography further: by comparing the model output with observations of snow depth, soil moisture and temperature at additional sites, though the study of the model's behaviour and uncertainty though testing a range of parameter values, and through gauging the impact of modelling microtopography on 865 modelled methane fluxes. We have shown that a simple two-tile representation of microtopography ( Figure 2) is able to improve modelled snow depths (Figure 7 landscape changes due to ground ice thaw and change in wetland area. The increase in methane fluxes was appreciable for 870 only the palsa-mire sites (+10 to 49%, Figure 13) due to the small difference in temperatures between tiles at the polygon sites ( Figure 11). We hypothesise that a similar effect on methane fluxes would be seen for other wetlands in the discontinuous permafrost region. This is because the model processes are general to all forms of microtopography, and because through this study we have found a weak dependence of soil temperatures on elevation differences greater than 0.5 m ( Figure C 1C, Appendix). To test this hypothesis would require data from additional sites, so this study underlines the need for more 875 microtopographically resolved data to be available. Snow depth was a key control on the modelled soil temperatures of each tile (Figure 12). Despite its simplicity, the snow scheme performed well at the modelled sites. While these sites are representative of widespread landforms, it remains to be seen if the scheme would perform as well in sites where the microtopography does not show such clear small-scale repetition. 880 Investigating the impact of linking the snow catch parameter to the vegetation height is a key next step and will be important for modelling snow-vegetation feedbacks. The impact of microtopography on modelled plant communities was not examined in this study as work needs to be undertaken to implement plant functional types in JULES which are more representative of arctic plants and their variety, and which have the correct response to saturated soils. The ability to model different plant functional types in microtopographically distinct areas would however be interesting with respect to peat formation and the 885 effect of different plant functional types on methane emissions (Cooper et al., 2017;Rupp et al., 2019). Combining this model with the ability to model peat formation as in Chadburn et al. (submitted) could lead to interesting dynamics as the tile elevations change, and would enable the future of features such as drained thermokarst lake basins to be modelled.
This study reiterates that making wetlands wet is a clear requirement for accurate representations of soil temperatures and 890 therefore methane emissions from high-latitude wetlands ( Figure 12). An explicitly modelled wetland 'tile', even if purely to provide more accurate soil temperatures for calculation of methane emissions, could therefore be a worthwhile step forward from models where gridcell soil temperatures are spatially homogeneous. Differences in the snowmelt, lateral drainage and runoff between tiles, as well as ponding on the low tile, were instrumental in modelling the effects of microtopography on hydrology. As such, it is difficult to suggest any of these processes that could be readily eliminated. However, several aspects 895 of the hydrology were identified which require consideration for future modelling. While much of our effort was directed at modelling unsaturated flows, it was found that these fluxes played only a small part in modelling these wetland environments.
Where lateral fluxes were appreciable, they were usually from a saturated layer to an unsaturated one. Modelling the lateral fluxes based on the position of the water table may therefore provide an acceptable basis for a simplified implementation of these flows, and one which may better model advective flows of heat. The amount of baseflow modelled from the low tile has 900 a large effect for the palsa mires. However, the current implementation of baseflow in JULES is not very suitable for a microtopographic model and also does not represent water influx, which is important for some wetlands. Consideration therefore needs to be given as to how the landscape-scale lateral flow out of and into wetlands could best be modelled. To allow for changes in the soil wetness, the methane flux calculation will need to be linked to the water table depth, which will also enable methanogen microbial dynamics to be better modelled. This will increase the sensitivity of the methane flux to the 905 water table depth that is already present due to the effect of soil moisture on soil temperatures. This sensitivity is part of the reason for the observed spatial variability of methane emissions from wetter areas. There is therefore a question of how upscaling the methane fluxes from those calculated for the low tile can be done in a representative way. Bechtold et al., (2019) suggest that for peatlands the TOPMODEL formulation of calculating the distribution of the water table should be replaced with a normal distribution based on the microtopographic variability, alongside other modifications to the hydrology. This 910 distribution could then be used alongside temperatures of tiles of different wetness to better represent the emissions from the whole of the mire. An alternative would be to explicitly represent the hydrological gradient using a larger number of tiles https://doi.org/10.5194/gmd-2021-285 Preprint. Discussion started: 26 October 2021 c Author(s) 2021. CC BY 4.0 License. and/or a representation of hillslopes (Nitzbon et al., 2020a;Swenson et al., 2019;Hazenberg et al., 2015). Both of these approaches could also provide a route to better modelling the landscape or catchment-scale baseflow. Currently the model does not include a treatment of surface waterbodies beyond being a simple surface store, for which an additional pond or lake 915 tile would be beneficial, as in Langer et al., (2016). Using a greater number of tiles may also be useful for the modelling of more complex features, for example by using three tiles to represent the centre, rim and trough in a polygonal landscape , or areas of land, wetland and lake (Schneider von Deimling et al., 2015). Indeed, Nitzbon et al. (2020) take this further, connecting multiple polygon units together to represent landscape-scale drainage and mesotopography.

920
Modelling microtopography provides a foundation for representations of thermokarst in Earth system models, and for answering key scientific questions about the impact of abrupt thaw processes on emissions and the future extent of wetlands.
Scaling to a pan-arctic model requires consideration of how model parameters such as tile areas and elevations can be assigned on this scale, and a representation of thermokarst requires a model of how these parameters can change over time. This study provides motivation for it being possible to assign parameters that are representative at scale, but also a measure of the 925 uncertainty as a result of estimating these parameters. In particular, we find that within the expected ranges at these sites, most parameters result in an associated uncertainty in July soil temperatures at 20 cm depth of under 0.8°, while the effect of the uncertainty in the high tile area and in the fractional drainage from the low tile result in an uncertainty on these temperatures of under 1.6° (Table 4). We also find a relatively weak dependence of summer soil temperatures on elevation difference within the expected ranges at each type of site. It would therefore be reasonable to associate a given elevation difference with 930 a particular landform. Previously, thaw models based on the thawing of excess ice and subsequent vertical elevation change have been used to simulate the future evolution of polygonal landscapes and palsas / peat plateaus (Langer et al., 2016;Aas et al., 2016Nitzbon et al., 2019Nitzbon et al., , 2020aCai et al., 2020;Martin et al., 2021). These approaches are particularly suited to polygonal landforms where the relative elevations of tiles determine the transformation between low centred and high centred polygons, affecting the hydrology of these landscapes (Nitzbon et al., , 2020a. However, for areas of peat plateau, and 935 for the expansion of thermokarst lakes, it may be more suitable to consider an area-based lateral thaw approach as this may require fewer tiles. We note that the modelled lateral heat conduction between tiles for Iškoras corresponds to melting 0.28 m 3 of ice per year per metre of palsa-mire border, which is of the same order of magnitude as the 0.13 m 3 yr -1 m -1 of volumetric loss of permafrost peat plateau edge found by Martin et al., (2021), indicating that such an approach may work. Recently, remote sensing and machine learning tools have enabled the mapping of the location and areas of different permafrost 940 landforms for large regions, and the mapping of how these are changing (Nitze et al., 2018). At the same time, climate envelope models have allowed for pan-arctic assessment of climatic suitability for palsas (Fewster et al., 2020) and ice wedge polygons, and which areas may be subject to change in the future. These data will enable microtopographic models to be initialised and provide a basis for model evaluation. The time is therefore ripe for explicit representations of microtopography to be set up and assessed on a pan-Arctic scale. 945

Code and data availability
Model code and Rose suite configuration files are freely available from the Met Office Science Repository Service ( https://code.metoffice.gov.uk/) upon registration and completion of a software licence. Details on accessing and running JULES configurations can be found in Wiltshire et al. (2020) Model output, processed observational data and plotting code is available at https://doi.org/10.5281/zenodo.5565163 (last access: 13 th October 2021). Sources for site observations can be found in Table A 2 (Appendix).

Appendices 1230
Appendix A -Site information and data sources