Interactive comment on “ Modeling the spatio-temporal variability in subsurface thermal regimes across a low-relief polygonal tundra landscape ”

J. Kumar et al. presented a pilot study using PFLOTRAN model to investigate the role of micro-topography in soil thermal dynamics of different types of ice-wedged tundra, which is important for further studying the responses of large-amount of frozen soil carbon to warming. Field measurements were provided for parameterization and validation. Therefore, I think the topic is important and the method is appropriate . The 3-D modeling is computing-intensive, it is hard, if not impossible, to be coupled in large-scale climate or terrestrial ecosystem models to investigate the effects of fine scale heterogeneity. Meanwhile, it is well-known that micro-topography of ice-wedged tundra ecosystem plays an important role in redistribution of surface water and vegetation growth. Therefore, the manuscript should focus on quantitatively assessing the


Introduction
Coastal Arctic landscapes -dominated by wetlands and patterned ground -cover approximately 5-10 % of Earth's land surface and play an important role in the hydrology, geomorphology, biogeochemistry and vegetation dynamics of the vast Arctic region.The low-gradient topography of the polygonal tundra characteristic of these landscapes is a complex mosaic of microtopographic features created by icewedge polygons.This microtopography leads to strong finescale variability in thermal and hydrological regimes of landscapes underlain by continuous permafrost.Permafrost landforms like drained lakes, low-centered polygons and highcentered polygons retard surface runoff after snowmelt, leading to increased surface water storage (in the form of lakes, ponds and wetlands) (Kane et al., 2003).Complex surface drainage patterns lead to heterogeneous soil moisture and substrate conditions, supporting a wide range of vegetation composition across the landscape.Arctic tundra soil pools are estimated to contain ∼ 190 Pg of carbon (Post et al., 1982), much of which is under risk of rapid release to the atmosphere in a warming climate.Hobbie et al. (2000) studied the controls over carbon storage and turnover in Arctic soils and found temperature, microtopography and vegetation composition to be the primary controls at regional scale.
Changes in the surface geomorphology which lead to the creation of ice-wedge polygons are induced by thermal disequilibrium and permafrost degradation.Lowland polygonal relief is dominated by low-centered polygons and highcentered polygons.Low-centered polygons are the most common polygonal landscape feature and are characteristic of poorly drained tundra.They consist of a raised rim with a wet central depression.Raised rims are the result of growing ice wedges that push material away from the center of the ice wedges to the sides (French, 2007).The standing water in the ice-wedge troughs leads to thermal erosion (i.e., accelerated thawing) along the rim.This preferential thaw may cause the ridge to collapse and form trenches surrounding the polygon center, inverting the relief to form high-centered polygons.High-centered polygons are well-drained with often dry centers, leading to low peat accumulation and deeper active layers.The microtopographic relief and associated heterogeneity in soil moisture support a diverse distribution of vegetation in the Arctic, with wet centers and troughs of low-centered polygons covered by mosses and sedges, while drier rims, and centers of high-centered polygons and rims of low-centered polygons are dominated by mosses, lichen and dwarf shrubs.These diverse land-cover types can also alter the surface energy balance and thermal properties through changes in albedo, surface roughness and evapotranspiration (Langer et al., 2011).
Large-scale climate and terrestrial ecosystem processes are represented in global-to regional-scale climate and ecosystem models.However, most of these models lack the representation of fine-scale heterogeneity in surface and subsurface processes at sub-grid scale that exercise significant control on the landscape-scale behavior.Representation of the fine-scale heterogeneity is important to model the nonlinear processes involved (Cresto Aleina et al., 2013).
Accurate characterization and modeling of subsurface thermal regimes in polygonal tundra is critical for our understanding of this sensitive system and our ability to predict its fate under climate change.In this study we developed approaches to (1) characterize the surface microtopography and subsurface structure of the polygonal tundra, (2) represent heterogeneous subsurface stratigraphy and hydraulic and thermal properties, (3) numerically model permafrost hydrology and (4) combine the above to simulate the permafrost thermal regime at field sites in a polygonal tundra region near Barrow, Alaska.

Study area
The study area is located within the Barrow Environmental Observatory (BEO), (Fig. 1) which lies 6 km east of Barrow, Alaska (71 • 18 N, 156 • 35 W), and is a field site of the US Department of Energy's Next Generation Ecosystem Experiments (NGEE) Arctic project.The BEO spans 32.21 km 2 of natural tundra, lakes and wetlands, and is reserved for scientific research.The landscape has low topographic relief, with elevations ranging from 0 to 7 m a.s.l.(above sea level) and low hydraulic gradient present across the region.Barrow has a polar maritime climate with mean annual air temperature of −12.0 and 3.3 • C during summer (June-August) (Liljedahl et al., 2011).The winter snowpack averages 20 to 40 cm, but snow accumulation is spatially variable due to variations in terrain roughness and drifting from strong easterly winds (Bockheim et al., 2001).Annual adjusted precipitation is 173 mm, with the majority of precipitation falling during summer months (Liljedahl et al., 2011).The polygonal tundra landscape is punctuated by thermokarst lakes and drained lake basins, with grass, moss and sedge as dominant vegetation types.Basins at Barrow are underlain by permafrost within 1 m of the surface and are classified as Gelisols, with an organic-rich surface layer underlain by a horizon of silt and clay to silt-loam textured mineral material and a frozen organic-rich mineral layer (Bockheim et al., 2001).The seasonal active layer thickness ranges between 30 and 70 cm at the BEO.
The various stages of geomorphologial and ecological change from low to transitional to high-centered polygons, lakes and drained lakes are all represented at the BEO.Fol- Relative elevation is the qualitative summary of topography in the region, while min/max/median are minimum, maximum and median elevations.lowing a "space for time" philosophy, NGEE Arctic intensive field sites at BEO were chosen across the landscape to observe and study polygonal landscapes at all stages of transition.Table 1 shows the characteristics of four sites (A, B, C, D) where our current study is focused.
A suite of observations is being collected at each of these intensive sites.Since 2012, meteorological data (including air temperature, summer precipitation, snow depth, relative humidity, wind speed and radiation) are being collected at all four sites (Hinzman et al., 2014b).Surface temperature data are being collected along a transect from the center of the polygon to ridge to trough.Each location on the transect consists of nine soil temperature sensors ranging in depth from 2 to 150 cm (Romanovsky and Cable, 2012).Figure C1 shows the observed time series of hourly air temperature and liquid precipitation at our four sites for the period of 1 October 2013-30 September 2014.Data for the 1-year period (1 October 2013-30 September 2014) were chosen for this study since it was the only complete year for which all the necessary observations were available.Figure 2a, d, g and j show high-resolution imagery of the sites where the observations have been collected.Figure 2a, d, g, and j show the boundary of the region around the intensive sites where we conducted the detailed modeling study presented here.

Methodology
To model the thermal regimes of the heterogeneous polygonal tundra ecosystem we developed approaches to (1) characterize the surface microtopography and subsurface structure of the polygonal tundra (2) represent heterogeneous subsurface stratigraphy and hydraulic and thermal properties and (3) numerically model permafrost hydrology.

Representation of landscape heterogeneity
Accurate representation of polygonal tundra in the model requires (1) identification of microtopographic features on the landscape and (2) characterization of soil stratigraphy and properties across the landscape.

Identification of polygonal features
The human eye can discern polygonal patterns and features in high-resolution satellite imagery with relative ease.However, automated recognition and delineation of polygonal features is challenging due to the variability in their spectral appearance, irregularity of polygon shape, dimension and orientation and lack of unique spectral signatures associated with the features (Skurikhin et al., 2013).Muster et al. (2012) investigated the subpixel heterogeneity in Landsat satellite imagery over ice-wedge polygonal tundra using a range of multi-scale data (field measurements and remote sensing), and concluded that resolutions of 4 m or less are necessary to map the fine-scale landscape elements of polygonal tundra.Skurikhin et al. (2013) used a combination of segmentation-and shape-based classification approaches using high-resolution WorldView-2 satellite imagery (60 cm resolution) to identify the landscape elements within the BEO.While they reported an overall accuracy of 95 %, their study region was limited to a 1000 × 1100 pixels subimage.The scalability of such a specialized algorithm based on high-resolution satellite imagery (of limited availability) is untested and difficult for application for landscape-scale studies like ours.
Thus, with landscape-scale application in mind, we employed a relatively simple and generic approach using a high-resolution digital elevation model (DEM) that exploits the relative difference in surface elevations that distinguish the polygonal features (center, ridges and troughs).Highresolution lidar data (25 cm resolution) were collected on 4 October 2005 by Tweedie (2010).The lidar data horizontal and vertical accuracy were approximately 30 and 15 cm respectively.Covering an approximately 2.5 km × 2.5 km area, the data set encompasses all of the NGEE Arctic intensive sampling sites (Fig. 1) where our study was focused.Using a high-resolution DEM created from this data set, elevation contours (10 cm interval) were developed to segment and classify the landscape in centers, ridges and troughs.
-Site D: Site D is relatively flat and is thus identified as flat-centered polygons, with the entire area within a narrow elevation range of 4.1-4.6 m.While polygonal features were evident in 0.25 m resolution aerial optical image (Fig. 2j), they were difficult to identify in the lidar DEM (Fig. 2k) due to the limitations of the vertical accuracy of lidar.Trough features in flat-centered polygons are not well pronounced.Thus the area was classified only as center (4.1-4.3 m) and tim (4.3-4.6 m) features.About 72 % of the area was classified as center while 28 % was classified as rim (Fig. 2l).
We did not apply any specialized rules to enforce any shape, dimension and/or patterns of the polygon features (center, ridge, trough), allowing us to scale our approach to the entire region where high-resolution DEMs were available.

Subsurface characterization
The structure and properties of subsurface soils are important factors controlling the pattern and variability of permafrost thermal processes in the tundra environment, and accurate characterization and representation of the heterogeneous subsurface properties is critical to understanding and modeling the subsurface thermal dynamics.However, the limited availability of soil properties in tundra environments and at our sites at the BEO presents a significant challenge.
During the period 31 July-3 August 2012 a field campaign was conducted by NGEE Arctic researchers to collect soil cores at one replicate polygon at each of the sites A, B, C and D at three microtopographic positions (center, ridge, trough) per polygon.Cores were collected using a hammer and a 5.08 cm diameter corer to collect one soil core per location to a depth of 30 cm.The soil horizons (moss, organic layer, mineral layer) were measured for each core to the nearest centimeter.A deep organic layer was found at several of the locations.However the total depth of the deep organic layer was not determined if it extended beyond the 30 cm core depth.Cores were collected for the purpose of biogeochemical analysis and thus no soil hydraulic or thermal properties were measured by the team.Figure 3 illustrates the subsurface soil horizons based on observations that we used in all our modeling studies presented here.In absence of colocated observations for soil hydraulic and thermal properties we derived the data for use in our studies from the published literature in tundra regions (Hinzman et al., 1991(Hinzman et al., , 1998) ) and a recent parameter calibration study conducted at one of our sites (Site C) (Atchley et al., 2015).Table D1 shows the soil hydraulic and thermal properties used in our modeling study.Collier and Kumar (2016) developed MeshMaker, a Pythonbased meshing framework to create high-resolution computational meshes for use in numerical simulation of permafrost thermal hydrologic processes at our polygonal tundra study sites.The meshing framework uses a high-resolution DEM and landscape classification (Fig. 2) to develop triangulated irregular networks (TINs).Nonuniform locally refined TINs adapt to the topographic complexity to create fine-resolution elements in areas with sharp changes in topography while creating coarser elements elsewhere, thus creating a highquality microtopography resolving mesh (Fig. 5).Variable resolution in vertical column was employed as described in Table 2. Data from Sects.3.1.1and 3.1.2were embedded within the generated meshes (Fig. 4) to represent the heterogeneity in the thermal hydrology models.By overlaying the TIN mesh with classified maps from Sect.3.1.1,microtopographic position (center/ridge/trough of a polygon) of each element in the mesh was identified.Polygon type and microtopographicspecific soil horizons data (Sect.3.1.2)were used to determine the soil horizons in the model mesh.While our data set was limited to a single replicate for each polygon type and location, significant spatial heterogeneity exists in reality.We assumed a variability of 10 % in soil horizon (moss, organic, mineral and deep organic soil) depths, and stochastically generated the soil horizon depths at each spatial location in the modeling domain.

Three-phase model for permafrost hydrology
In this study, we will use the open-source code PFLOTRAN to model the flow of mass and energy in the subsurface.PFLOTRAN (Hammond et al., 2016(Hammond et al., , 2014) ) is a state-of-theart, massively parallel subsurface flow and reactive transport code.PFLOTRAN solves a system of generally nonlinear partial differential equations (PDEs) describing multiphase, multicomponent and multiscale reactive flow and transport in porous materials.The PDEs are spatially discretized using a finite volume technique, and the backward Euler scheme is used for implicit time discretization.One system of PDEs that PFLOTRAN implements is a three-phase, thermal hydrology model (the TH process model in PFLOTRAN parlance), which describes a balance of mass and energy in which the liquid pressure P and the bulk temperature T are the primary unknown variables.In Eqs. ( 1) and ( 2), ϕ refers to porosity, s to percent saturation, η to molar density, U to internal energy, ρ to mass density, c to specific heat and H to enthalpy.The subscripts { , g, i} refer to the liquid, gas and ice phases of water, respectively, and the subscript "s" to the soil matrix.The Darcy velocity is given by where k denotes intrinsic permeability, k r relative permeability, µ viscosity, g unsigned gravity and z the vertical component of the position vector x.The effective thermal conductivity is expressed as where κ are the thermal conductivities of each pure phase and K represents the Kersten number of the frozen and unfrozen phase: where = 1 × 10 −6 and α i and α are parameters of the assumed power law. Figure 6 shows the modeled dependence of effective thermal conductivity (κ eff ) on the fraction of water present in liquid (s )), ice (s i ) and gas (s g ) phases.
The variables Q M and Q E represent generic mass and energy sources and sinks.We emphasize that the saturations, densities and internal energies are all nonlinear functions of the liquid pressure and temperature and include the latent heat of fusion associated with change of phase.We also note that PFLOTRAN implements several choices of constitutive models for relating the saturations to the liquid pressure and bulk temperature.While only a brief overview of the numerical formulation has been provided here, we would refer to Painter (2011), Painter and Karra (2014) and Karra et al. (2014) for detailed discussion of the formulation.Key parameters for the model relevant for current study are described in Table D1.

Initial and boundary conditions
Three-dimensional subsurface models for each of the four sites were initialized by freezing the entire modeling domain at a temperature of −1.0 • C. The models were spun up to a thermal periodic steady state using a time series of mean daily temperatures applied to the top of the domain (ground surface).Spin-up simulations were conducted for a period of 10 years by cycling annual time series of forcing.Spin-up simulations were continued until a periodic steady state was achieved (i.e., close to zero interannual variability in annual thermal regime).Spin-up duration of 10 years was used at all the sites and was determined to be sufficient.We conducted a series of initialization simulations by varying initial temperatures at start of spin up and found them to not have any significant impact on the final periodic steady state, besides simulation period required to reach that steady state.Mean daily near-surface temperature time series for the period 1 October 2013-30 September 2014 were derived from hourly in situ temperatures from sensors located at 2 cm depths.At all four sites, using sensors installed at center, ridge and troughs, three different time series were prepared.Using the classification (center, ridge, trough) embedded in the model (Sect.3.1.1),these microtopography-specific temperature time series were applied in a spatially heterogeneous, microtopography-aware fashion to simulate the complex thermal hydrologic regimes in permafrost soils at the BEO.A no-flow boundary condition was applied to the sides of the domain, while the deep bottom boundary was held at a constant temperature of −10 • C, based on the temperature from the West Dock site (Fig. 3 of Romanovsky et al., 2010), which is located at a comparable latitude (70.4 • N) to the BEO (71.29 • N).It would be important to note that due to lack of observations of drainage pattern at the site, we opted to use a no-flow boundary condition; however, surface runoff in and out of the region occurs in reality.Thus the soil moisture states were not very well constrained in the model, and this also has consequences for the simulated thermal regime.Surface processes (such as vegetation cover and snow) play an important role in regulating the thermal regimes of permafrost soils.While the surface processes are not represented in our model, use of near-ground surface (2 cm depth) temperature as the boundary condition for the simulation allows us to isolate (though not completely) the effect of surface processes.

Simulation of permafrost thermal regimes
After the models were spun up to periodic steady-state condition, the simulation was continued for another year, and outputs were used for validation and analysis.Soil temperature observations from the thermal sensors at 2 cm depths at the sites for the period 1 October 2013-30 September 2014 were used to drive the time-dependent (Dirichlet) boundary condition at the top (ground surface) of the model.In addition, mean daily time series (1 October 2013-30 September 2014) of liquid (summertime) precipitation was also applied as moisture input to the model.Groundwater infiltration was considered to be zero if ground surface temperature was below freezing or if the domain was fully saturated.While soil moisture plays an important role, the focus of this study was on the thermal regime, and thus all results and discussions presented are focused on soil temperature.Section E3 shows the spatial pattern of maximum water table depth across the study region and exhibits strong correlation with microtopography.Soil temperature data from all the sensors (16 sensors at depths from 5 to 150 cm) were used to evaluate the accuracy of the models.While we have selected sensors at four depths of 5, 10, 50 and 150 cm for discussion, the results at all depths can be found in Sect.E1.
At all the sites, simulated soil temperatures in the topmost soil layer (5 cm thick) compared well with the observed nearsurface temperature at 5 cm depth (Figs.7a-10a).Simulated soil temperatures in deeper soils show warm bias increasing with depth as compared to the observed temperatures at the sensors (coefficient of determination R 2 of 0.93-0.99)(Tables 3-6).The model matched the observations with a root mean square error (RMSE) of 0.60-0.99• C near the surface, with an increasing error in deeper soils.While the modeled temperature bias was in the range of −0.30-0.10• C near the surface, a warm bias of up to 1.0-1.8• C was found in deep soils.Warm bias in the model was more pronounced during summer season, compared to winter season.
Figures 7-10 show the comparison between the simulated and observed soil temperatures at several select depths (5, 10, 50 and 150 cm from surface) at sites A, B, C and D respectively.Simulated temperatures across all the sites matched well with the observed temperatures at shallow depths, but showed a deviation towards temperatures that are warmer than observed in deep soils.A number of factors may be contributing to this bias in simulations, including applied boundary conditions (Sect.3.3).We believe that insufficient characterization and parameterization of heterogeneous properties due to limited data availability is one of the key reasons for this bias.The model was evolved to a periodic thermal steady state through a spin-up process, and soil temperatures in vertical profile are strongly dependent on the soil properties.Forced with surface temperature boundary conditions, while the soil temperature in the model after the spin-up stage The Cryosphere, 10, 2241-2274, 2016 was close to that observed, a warm bias was observed in the deeper soils.That warm bias was carried over to the final stage of the simulation, resulting in bias in the simulated soil temperatures and thaw depths reported here.With soil cores collected at the sites limited to the top 30 cm of the soil, our understanding of structure and physical and thermal properties of deeper soils is limited.For example, while we know that the presence of ground ice (like ice wedges, segregated ice, ice lens etc.) is common in the subsurface of Arctic tundra, their representation in the model is completely missing.A lack of representation of these cryostructures is potentially one of the reasons for warmer soils in our simulations.While PFLOTRAN has the ability to capture and model such cryostructures (via heterogeneous subsurface structure and properties but not their formation and evolution), we lack any quantitative data to characterize them for representation in the model.Ongoing efforts under the NGEE Arctic project by Kneafsey and Ulrich (2016)  Heat flow in the permafrost soils, which are frozen for a significant part of the year, occurs primarily due to conduction.Thermal conductivity of the soil is sensitive to the temperatures and thus the fraction of water present in liquid vs. ice phase (Eq.4, Fig. 6).Effective thermal conductivity (κ eff ) of the soil is higher during the winter months when almost the entire soil domain is in a frozen state (thus high ice saturation -s i -and low liquid saturation -s ), compared to summer months when active layer is in a thawed state.Higher conductivity of the soil, and thus higher conductive heat flows during the winter season in a heterogeneous soil domain, leads to higher temperature variability in model simulations of the permafrost thermal regimes at our sites.
During the summer season, advective heat flow processes occur within the thawed soil layers.The study sites also receive liquid precipitation and thus infiltration during the summer season (Fig. C1), which leads to vertical as well as horizontal flows in the thawed soil.For example, at Site A horizontal velocities are close to zero during the early summer (Fig. 12a) when soil temperatures are close to freezing.Af-    ter the ground has thawed, liquid precipitation events during summer (Fig. 12b and c) lead to significant lateral flows.High-elevation rim regions drain to the center and trough of the low-centered polygons.
Large summer bias in the model is partly due to the illconstrained flow boundary conditions and the lack of explicit representation of surface flow processes.As ground thaws with the rise of summer temperature, significant lateral flow occurs, even in the absence of rainfall events (Fig. 12d).No flow boundary condition applied in our model prevents any runoff out of the modeling domain, leading to increased soil moisture (Fig. E4) in low-elevation areas where the higher warm bias in the model is observed .
A reference 1-D simulation was conducted at Site A (described in Sect.E4) to isolate and understand the role of lateral flows.Simulations show that ignoring lateral flows can lead to a temperature bias of −0.55-0.85• C across the site depending on the microtopography (Fig. E5).
In a permafrost environment, the active layer is the top layer of soil that thaws during the summer and freezes again during the autumn.While PFLOTRAN solves for soil temperature as the primary variable, we derived active layer     depth (or thaw depth) as the sum of the thickness of soil layers above freezing temperature (0 • C).Figures 13-16 show the temporal dynamics of thaw depth during the period of simulation and spatial pattern of maximum thaw depth across the region.A wide spatial variability in the thaw depth was observed in the simulations, which were strongly correlated to the microtopography.The variability is primarily derived by the microtopography and the subsurface heterogeneity.The warm bias in soil temperature in the model translates to a bias towards deeper thaw depths (0.60-1.0 m) as compared to the observations (0.30-0.74 m) at the site (Sect.F and Fig. F1).
One-dimensional reference simulations (Sect.E4) demonstrate the importance of 3-D representation of thermal and hydrologic processes to capture the effect of microtopography.While the lateral flow may not significantly effect the average thaw depths across a large region (Fig. E7a), they lead to microtopography-dependent differences in thaw depths (Fig. E7b-d) throughout the thaw season, which may have an implication for other ecological processes like biogeochemistry and vegetation dynamics.

Understanding the thermal regimes of polygonal tundra
Microtopography of the polygonal tundra exerts critical controls on the flow of water and energy at local to regional scales, which further influences the ecological and biogeo-chemical processes on the landscape.Surface processes (not modeled in our study) like vegetation and snow cover also play a critical role in regulating the subsurface thermal regimes through thermal insulation effects.In our modeling approach we represented the microtopographic features center, rim and trough across four low to transitional to highcentered polygons.PFLOTRAN successfully simulated the overall pattern of thermal regimes in center, rim and trough across the four sites A, B, C and D, though a significant warm bias is observed in deeper soils (Table 1).
-Site A: Site A is located in a poorly drained region dominated by low-centered polygons with lowelevation centers, raised rims and troughs.Center areas are warmer than rim and trough areas, while rims are coldest (Fig. 7).Centers in low-centered polygon are often inundated and relatively wet (Fig. E4) most of the year and support vegetation (mosses and sedges).
Low-elevation centers also receive higher snow cover.Vegetation and snow cover provide thermal insulation to the ground, keeping the center region warmer compared to rim and trough.Dry rims (Fig. E4) with low vegetation cover and low snow accumulation (Fig. G1) are most exposed to the winter temperatures, and are thus the coldest.
Three-dimensional representation of the processes is important to simulate these microtopographydriven processes.One-dimensional reference simulation (Sect.E4) shows that the lack of full 3-D representation can lead to topography-driven bias in simulated thermal regimes (Fig. E7).
-Site B: Site B is dominated by well-drained highcentered polygons with relatively dry elevated centers and deep troughs (Fig. E4).Vegetation in high-centered polygons are dominated by lichens, moss and dwarf shrubs.In contrast to the low-centered polygons, centers and dry-tundra graminoids, which have low vegetation and snow cover (Fig. G1), are most exposed to the changes in air temperatures, and are thus colder than rim and trough, which show relatively warmer soil temperature regimes during the winter (Fig. 8).
-Site C: Site C is located in an area of geomorphological transition from low-to high-centered polygons, characterized as flat-centered polygons.They consist of shallow flat centers, deep troughs and raised rim regions.Soils in deep troughs are thermally insulated by higher snow cover (Fig. G1) and thus show warmer soil temperature regimes compared to centers and rims (Fig. 9).Center and rim regions show similar thermal regimes with centers being slightly warmer due to higher vegetation and snow cover.
-Site D: Site D is characterized by low-centered polygons with no pronounced rims.Site D is wettest among the four study areas, with low-elevation center areas that remain inundated for most of the summer season (Fig. E4).While snow accumulation was fairly uniform (Fig. G1) across the flat region, vegetation cover plays an important role.Wet centers support rich vegetation, leading to warmer soil temperatures as compared to the trough regions (Fig. 10).
5 Model uncertainties and limitations

Why not to calibrate
Accurate simulation of permafrost thermal regimes requires the mechanistic representation of thermal hydrologic processes in the model.However, equally important is the accurate representation of subsurface structure and soil properties, model parameters and initial and boundary conditions.Given the lack of co-located observations for the soil properties, in this study we used soil thermal and hydraulic properties data from a different tundra site based on Hinzman et al. (1998).While PFLOTRAN was able to simulate the thermal hydrologic processes and match the soil temperature observations at the sites across a range of polygonal landscape and microtopography features fairly well, simulated soil temperatures show deviations from the observed temperatures at times.Simulated temperatures show warm bias in deep soils where data for soil characterization and properties are almost completely missing.Parameter calibration is a popular technique that has been widely used in hydrologic modeling to determine model parameters and properties to optimize the model fit to target observations.The high-resolution 3-D PFLOTRAN thermal hydrology model used in this study includes many degrees of freedom and parameters, which combined with the complex nonlinearity of hydrologic processes poses a complex high dimensional optimization problem.While a wide range of calibration approaches (Heuvelmans et al., 2006;Madsen, 2000;Shafii and De Smedt, 2009;Singh and Minsker, 2008) are available to determine optimal model parameters to fit the observed data (soil temperatures in this study), we face the problem of nonuniqueness (equifinality).A diverse set of possible parameter values can lead to similar model performance (Beven and Freer, 2001).The issue of non-uniqueness is especially pronounced in tundra ecosystem due to poor availability of data and thus poor bounds on parameters, which leads to a high degree of uncertainty in the models with or without calibration.
While systematic calibration can help identify effective parameters for the model, transfer of parameters across models and modeling domains is difficult (Bárdossy, 2007).At data-limited study sites like ours, while a ill-constrained calibration may compensate for lack of data, it does not improve our understanding of the system.In this study we choose not to calibrate the model parameters to achieve better fit with The Cryosphere, 10, 2241-2274, 2016 the observations, instead using the uncalibrated results to diagnose the potential model deficiencies and identify the characterization and data needed to better represent the real world in our simulations.
When properly constrained using observed data, the model calibration is a powerful tool that can provide improved physical representation of processes and parameters and new insights.Atchley et al. (2015) and Harp et al. (2016) have successfully demonstrated the use of such techniques for three-phase hydrology models at one of our sites (site C) and would inform our future studies.

Identifying model and data gaps
While the agreement between modeled and observed simulated soil temperature demonstrates the ability of the model to simulate the thermal hydrologic processes in the polygonal tundra, disagreements help us identify the existing gaps in data and model.
Modeling results highlight the need for co-located measurements of soil thermal and hydraulic properties for accurate modeling of hydrologic processes.While most soil core observations, including those used in our study, are focused on the shallow active layer, characterization of deeper permafrost soils is essential for understanding the thermal regimes and potential changes expected under warming climate.Modeling soil temperatures, beyond the high-level estimation of thaw depth (or active layer thickness), is important to understand the thermal regime of permafrost soil and its behavior under warming conditions.For example, during winter seasons even when the soils are completely frozen, variability in soil temperatures (Figs.7-10) exists and may impact carbon fluxes from the system (Zona et al., 2016).Soil temperatures that are warmer than observed in deep soils in our thermal periodic steady-state solutions are due to inaccurate soil characterization and poorly bounded boundary conditions at the bottom of the modeling domain.Heat flux observations in deep permafrost, while hard to measure, would help provide accurate bounds for the thermal hydrology model.In addition to rainfall events (which were captured in our models), surface drainage processes provide inputs to the ground water system (not captured in our models).Surface drainage observations in the local catchments (not available to us) are needed to appropriately model and constrain this process.

Summary and conclusions
Low-relief polygonal tundra ecosystems consist of microtopographic features that controls the local-scale hydrology.The water and energy flow patterns on the landscape in turn regulate biogeochemical processes and vegetation dynamics.The objective of this study was to develop an end-toend modeling approach for landscape-scale modeling of per-mafrost thermal hydrology of real-world polygonal tundra sites, improving our ability to model and to understand the patterns of thermal regimes.Using the best available data for our study sites we developed techniques to characterize the polygonal microtopographic features and represent the heterogeneous soil hydraulic and thermal properties at the sites.These data sets were embedded within topography-following high-resolution meshes to simulate the thermal hydrologic processes in PFLOTRAN thermal hydrology model.We employed detailed surface meteorology and subsurface soil temperature observations from the site to simulate and analyze the thermal regimes at four representative sites in polygonal tundra ecosystem at the Barrow Environmental Observatory.
Our modeling-based study reveals the role of microtopographic features in regulating the permafrost thermal dynamics across heterogeneous polygonal tundra landscape.Simulation results at four sites across the polygonal tundra landscape demonstrate the viability of the developed mechanistic approach to model the thermal regimes.Thermal regimes of center, rim and trough features of polygonal tundra exhibit distinct patterns in low to transitional to high-centered polygon landscapes, which are governed by the microtopography, surface and subsurface hydrology and surface processes (like air temperature, snow cover and vegetation).Our reference 1-D model simulation demonstrates that while 1-D approaches may be sufficient to understand mean thermal behavior of a large landscape, detailed 3-D process representation is required to capture microtopography controls at high resolution.Our PFLOTRAN-based modeling approach was able to simulate these overall patterns at four study sites.Comparing the simulated soil temperatures to available observations, the model demonstrates the ability to simulate the soil temperature at shallow depths, but deviations from observations in deep soils highlight the need for better soil characterization using deep cores in these ecosystems.Our study also highlights the need for co-located observations for accurate modeling and understanding of the tundra landscape.Model disagreements with the observations in this study may partially be due to use of soil properties from literature in the absence of site-based measurements.Under the NGEE Arctic project, we are working with field scientists for improved co-located measurements.The present study was limited to a single year for which we had all the necessary data for model forcings and validation available, and thus was not able to address the role of interannual variability.We plan to address this important problem as more data from our sites become available.While we have not addressed all the deficiencies in model process representation and parameterization identified and reported here in this study, we believe we have developed and presented a process-rich modeling framework as a first critical step that would enable such studies.The modeling approach developed in this study -when combined with necessary co-located observations -promises to allow accurate modeling of permafrost thermal hydrology and provide a framework for identifying and guiding the future observa-tions required for improved modeling and understanding of the polygonal tundra ecosystem.
In a warming world, wet low-centered polygon landscapes are expected to go through geomorphological change to drier high-centered polygon landscapes.The modeling approach developed in this study is a first step towards enabling future investigations of the impact of thermal hydrologic changes in these landscapes under projected climate scenarios.The model developed here does not have the ability to simulate the dynamic changes in microtopography expected due to ice-wedge degradation (Liljedahl et al., 2016).While beyond the scope of the current study, ongoing developments in biogeochemical modeling within PFLOTRAN (Tang et al., 2015) in combination with our thermal hydrology model developments will also allow modeling of the terrestrial carbon cycle in this sensitive landscape under future warming scenarios.
While the knowledge gained by developing and evaluating fine-scale 3-D simulations is valuable from the perspective of increased understanding of complex process interactions, the explicit long-term goal of the NGEE Arctic project is to improve predictions of Arctic ecosystem processes at scales relevant to coupled climate and Earth system simulation.One element of our strategy to migrate knowledge across scales is to improve the grid and sub-grid representations in the land model component of our Earth system models to capture observed modes of variability in physical, biological and biogeochemical processes.For example, our new top-level grid topology for global-scale land modeling follows watershed boundaries instead of the typical and arbitrary rectangular grid cell arrangement (Tesfa et al., 2014).Sub-grid schemes are being developed that represent topographic variation within basins, and our goal is to apply these methods in the microtopographic setting of polygonal tundra to capture the variation in thermal, hydrologic and biogeochemical regimes, and interactions with vegetation communities.The current study is one step toward identifying the relevant modes of variation among diverse landforms in the polygonal tundra region.Another element of our scaling strategy is to use, to the full extent possible, a common set of modeling tools to construct simulations at various spatial scales.Even though many processes that can be represented explicitly at the finest scales (such as lateral flows of energy and water) must be parameterized for efficiency in a larger scale simulation, having a common underlying set of equations helps to reduce unintentional loss of information across scales due, for example, to aggregation and disaggregation operators.

Data availability
Model source codes, input, forcings and output data sets from the study are publicly available at doi:10.5440/1184018(Kumar et al., 2016).To evaluate the benefits of 3-D representation of thermal and hydrologic processes over 1-D representation, a reference 1-D simulation was conducted at Site A. To allow for a fair comparison of the 3-D vs. 1-D models, we implemented a vertical-flow-only mode in PFLOTRAN that allows all lateral flows to be turned off.The mode was applied at Site A with exactly the same mesh, soil properties, initial conditions and forcings as the 3-D simulation at Site A (Sect.4.1), effectively converting the 3-D modeling domain into 2305 1-D columns.

Appendix D: Soil hydraulic and thermal properties
A lack of horizontal flows in 1-D simulations shows an impact on the soil temperatures during the summer season, when lateral flows are observed in the 3-D model (Fig. 12).Differences in temperatures are strongly correlated with the topography of the domain, with high-elevation areas experiencing warm bias (water and energy not leaving due to lack of horizontal flow), while low-elevation areas experience cold bias (water and energy not entering the cell due to lack of horizontal flow) (Fig. E5).A temperature difference ranging −0.55-0.85• C was observed across Site A during the summer season.
A lack of horizontal flow also affects the simulated estimates of thaw depth across the study region.A difference in maximum thaw depth of −10 to +20 cm was observed across the entire study region during the summer season (Fig. E6).These differences are primarily due to the lack of lateral flows in 1-D models and demonstrate the strong effect of microtopography.High-elevation areas (like rims) do not lose water and energy due to lack of lateral flows, and show increased thaw depths, while low-elevation areas (like polygons' centers) show reduced thaw depths.

Figure 1 .
Figure 1.NGEE Arctic field sites at BEO (lidar elevations are in units of meters above mean sea level).

Figure 2 .
Figure 2. Elevation contour-based classification of the study areas at Sites A, B, C and D. In subfigures (c, f, i, l) the colors reflect polygon type (red: center; green: ridge; blue: trough).

Figure 3 .
Figure 3. Subsurface soil horizons across microtopographic positions at Sites A, B, C and D based on field observations.Models at study sites were parameterized using these data.

Figure 4 .
Figure 4. Workflow for heterogeneous soil parameter assignment in the computational mesh.

Figure 5 .
Figure 5. Microtopography resolving unstructured meshes.The coloring reflects the ground surface elevation.Variable resolution mesh adapts to complex topography and uses finer resolution in areas of strong relief and coarser resolution in simpler terrain.

Figure 6 .
Figure 6.Model dependence of effective thermal conductivity on liquid-(s ), ice-(s i ) and gas-phase (s g ) fraction of water.

Figure 7 .
Figure 7. Simulated vs. observed soil temperatures at Site A at depths 5, 10, 50 and 150 cm show bias in the model in deeper soil and during the warm summer season.Dotted lines represent observed data at center (in red), rim (green) and trough (blue) locations, while shaded curves show mean ± standard deviation of simulated daily soil temperatures across the domain.
and Dafflon et al. (2016) using X-ray computed tomography (CT) scanner technology on ice cores from BEO can potentially provide detailed 3-D soil structure and density information and to help address this missing piece.Spatial variability in soil temperatures was observed in the simulations (Figs.7-10), arising in part due to 3-D heat flow and heterogeneous subsurface structure and soil properties represented in the model.Simulated soil temperatures also show a seasonal pattern of spatial variability with high vari-ability during the cold winter season and lower spatial variability during summer.Figure 11 shows the time series of spatial variability (standard deviation) in soil temperature at Site A during the simulation period, showing strong seasonality, the magnitude of which is reduced at deeper soils.Similar patterns of variability were observed at Site B and C, depicted in Figs.E1-E3).

Figure 8 .
Figure 8. Simulated vs. observed soil temperatures at Site B at depths 5, 10, 50 and 150 cm show bias in the model in deeper soil and during the warm summer season.Dotted lines represent observed data at center (in red), rim (green) and trough (blue) locations, while shaded curves show mean ± standard deviation of simulated daily soil temperatures across the domain.

Figure 9 .
Figure 9. Simulated vs. observed soil temperatures at Site C at depths 5, 10, 50 and 150 cm show bias in the model in deeper soil and during the warm summer season.Dotted lines represent observed data at center (in red), rim (green) and trough (blue) locations, while shaded curves show mean ± standard deviation of simulated daily soil temperatures across the domain.

Figure 10 .
Figure 10.Simulated vs. observed soil temperatures at Site D at depths 5, 10, 50 and 150 cm show bias in the model in deeper soil and during the warm summer season.Dotted lines represent observed data at center (in red), rim (green) and trough (blue) locations, while shaded curves show mean ± standard deviation of simulated daily soil temperatures across the domain.

Figure 11 .Figure 12 .
Figure 11.Spatiotemporal variability (standard deviation) in simulated soil temperature at Site A at depths 5, 10, 50 and 150 cm.High variability in soil temperatures was observed across the region during the year, especially cold winters with centers showing highest variability.Red, green and blue lines represent center, rim and trough respectively, and the black line represents the standard deviation across Site A.
Figure 14.Temporal (left panels) patterns of simulated thaw depths (defined as thickness of soil layers above 0 • C) and spatial pattern of maximum thaw depth (right panels) at Site B. The bold line represents the mean thaw depth across the site, while the shaded curve represents standard deviation.

Figure E1 .
Figure E1.Spatiotemporal variability (standard deviation) in simulated soil temperature at Site B. Red, green and blue lines represent center, rim and trough respectively, and the black line represents the standard deviation across Site B.

Figure E2 .
Figure E2.Spatiotemporal variability (standard deviation) in simulated soil temperature at Site C. Red, green and blue lines represent center, rim and trough respectively, and the black line represents the standard deviation across Site C.

Figure E3 .
Figure E3.Spatiotemporal variability (standard deviation) in simulated soil temperature at Site D. Red, green and blue lines represent center, rim and trough respectively, and the black line represents the standard deviation across Site D.

Figure
FigureE4shows the spatial distribution of maximum water table elevation during the simulation period at sites A, B, C and D. Ground surface is at 50 m elevation at all the sites.

Figure E4 .
Figure E4.Elevation of maximum water table across the study region.Ground surface is at 50 m elevation.The bold line represents the mean thaw depth across the site, while the shaded curve represents standard deviation.

Figure E5 .
Figure E5.Difference in mean ground surface temperature during the early summer season in 1-D vs. 3-D models.Compared to 3-D models (Fig. 12), temperature differences observed in 1-D models are due to the lack of horizontal flows.Positive differences indicate warmer conditions in 1-D models compared to 3-D models, while negative temperatures indicate cold bias.

Figure E6 .
Figure E6.Difference in maximum thaw depth between 1-D and 3-D simulations.A positive bias indicates higher thaw depth in 1-D simulations compared to 3-D models, while negative values indicate lower thaw depths.

Figure E7 .
Figure E7.Mean difference in thaw depth between 1-D models compared to 3-D models.A positive bias indicates higher thaw depth in 1-D simulations compared to 3-D models, while negative values indicate lower thaw depths.Daily precipitation patterns (in millimeters per day) are shown by blue bars.

Table 1 .
Areas A, B, C, D polygonal features and environmental characteristics.

Table 2 .
Mesh resolution in vertical column.Finer resolution was used in the layers closer to the surface, while a coarser resolution was used in deeper soils.

Table 3 .
Model performance statistics at Site A compared to observed soil temperatures (for annual, winter and summer periods).The model shows higher bias during the summer season compared to winter.
RMSE is the root mean squared error and R 2 is the coefficient of determination; a negative bias indicates cold bias in the model, while a positive bias indicates a warm bias.

Table 4 .
Model performance statistics at Site B compared to observed soil temperatures (for annual, winter and summer periods).The model shows higher bias during the summer season compared to winter.
RMSE is the root mean squared error and R 2 is the coefficient of determination; a negative bias indicates cold bias in the model, while a positive bias indicates a warm bias.

Table 5 .
Model performance statistics at Site C compared to observed soil temperatures (for annual, winter and summer periods).The model shows higher bias during the summer season compared to winter.
RMSE is the root mean squared error and R 2 is the coefficient of determination; a negative bias indicates cold bias in the model, while a positive bias indicates a warm bias.

Table 6 .
Model performance statistics at Site D compared to observed soil temperatures (for annual, winter and summer periods).
RMSE is the root mean squared error and R 2 is the coefficient of determination; a negative bias indicates cold bias in the model, while a positive bias indicates a warm bias.

Table D1 .
Soil hydraulic and thermal properties used in the models.

Table E1 .
Model performance statistics at Site A compared to observed soil temperatures.
RMSE is the root mean squared error and R 2 is the coefficient of determination; a negative bias indicates cold bias in the model, while a positive bias indicates a warm bias.

Table E2 .
Model performance statistics at Site B compared to observed soil temperatures.
RMSE is the root mean squared error and R 2 is the coefficient of determination; a negative bias indicates cold bias in the model, while a positive bias indicates a warm bias.

Table E3 .
Model performance statistics at Site C compared to observed soil temperatures.is the root mean squared error and R 2 is the coefficient of determination; a negative bias indicates cold bias in the model, while a positive bias indicates a warm bias. RMSE

Table E4 .
Model performance statistics at Site D compared to observed soil temperatures.