Impacts of microtopographic snow-redistribution and lateral subsurface processes 1 on hydrologic and thermal states in an Arctic polygonal ground ecosystem 2 3

Microtopographic features, such as polygonal ground, are characteristic sources of landscape heterogeneity in the Alaskan Arctic coastal plain. Here, we analyze the effects of snow redistribution (SR) and lateral subsurface processes on hydrologic and thermal states at a polygonal tundra site near Barrow, Alaska. We extended the land model integrated in the ACME Earth System Model (ESM) to redistribute incoming snow by accounting for microtopography and incorporated subsurface lateral transport of water and energy (ALMv0-3D). Three 10-years long simulations were performed for a transect across polygonal tundra landscape at the Barrow Environmental Observatory in Alaska to isolate the impact of SR and subsurface process representation. When SR was included, model results show a better agreement (higher R 2 with lower bias and RMSE) for the observed differences in snow depth between polygonal rims and centers. The model was also able to accurately reproduce observed soil temperature vertical profiles in the polygon rims and centers (overall bias, RMSE, and R 2 of 0.59 °C, 1.82 °C, and 0.99, respectively). The spatial heterogeneity of snow depth during the winter due to SR generated surface soil temperature heterogeneity that propagated in depth and time and led to ~10 cm shallower and ~5 cm deeper maximum annual thaw depths under the polygon rims and centers, respectively. Additionally, SR led to spatial heterogeneity in surface energy fluxes and soil moisture during the summer. Excluding lateral subsurface hydrologic and thermal processes led to small effects on mean states but an overestimation of spatial variability in soil moisture and soil temperature as subsurface liquid pressure and thermal gradients were artificially prevented from spatially dissipating over time. The effect of lateral subsurface processes on active layer depths was modest with mean absolute difference of ~3 cm. Our integration of three-dimensional subsurface hydrologic and thermal subsurface dynamics in the ACME land model will facilitate a wide range of analyses heretofore impossible in an ESM context.

Several recent modeling studies have predicted a positive carbon-climate feedback at the global scale (Cox et al., 2000;Dufresne et al., 2002;Friedlingstein et al., 2001;Fung et al., 2005;Govindasamy et al., 2011;Jiang et al., 2011;Jones et al., 2003;Koven et al., 2015;Matthews et al., 2007b;Matthews et al., 2005;Sitch et al., 2008;Thompson et al., 2004;Zeng et al., 2004), although the strength of this predicted feedback at the year 2100 was shown to have a large variability across models (Friedlingstein et al., 2006).In contrast to the ocean carbon cycle, the terrestrial carbon cycle is expected to be a more dominant factor in the global carbon-climate feedback over the next century (Matthews et al., 2007a;Randerson et al., 2015).
Changes in Arctic ecosystem net ecosystem productivity (NEP, defined as the difference between net primary production (NPP) and heterotrophic respiration (Rh)) will be determined by the magnitude and direction of changes in NPP and Rh.Warming experiments in the Arctic have found increases and decreases of plant growth in response to higher temperatures (Barber et al., 2000;Chapin et al., 1995;Cornelissen et al., 2001;Hobbie and Chapin, 1998;Hollister et al., 2005;Van Wijk et al., 2004;Walker et al., 2006;Wilmking et al., 2004).Arctic ecosystems are limited in nitrogen availability (Schimel et al., 1996;Shaver and Chapin III, 1986) and higher mineralization rates under warmer climate (Hobbie, 1996) could lead to higher CO2 fixation by plants (Shaver and Chapin, 1991).
Additionally, a longer growing season is expected to result in a negative carbon-climate feedback by increasing NPP (Euskirchen et al., 2006).On the other hand, microbial decomposition of previously frozen soil organic matter under a warmer climate is expected to strengthen the carbon-climate feedback (Davidson and Janssens, 2006;Mack et al., 2004;Oechel et al., 1993;Tarnocai et al., 2009).Snow, which covers the Arctic ecosystem for 8-10 months each year (Callaghan et al., 2011b), is a critical factor influencing hydrologic and ecologic interactions (Jones, 1999).Snowpack modifies surface energy balances (via high reflectivity), soil thermal regimes (due to low thermal conductivity), and hydrologic cycles (because of melt water).
Several studies have shown that warm soil temperatures under snowpack support the emission of greenhouse gases from belowground respiration (Grogan and Chapin Iii, 1999;Sullivan, 2010) and nitrogen mineralization (Borner et al., 2008;Schimel et al., 2004) during winter.Additionally, decreases in snow cover duration have been shown to increase net ecosystem CO2 uptake (Galen and Stanton, 1995;Groendahl et al., 2007).Recent snow manipulation experiments in the Arctic have provided evidence of the importance of snow in the expected responses of Arctic ecosystems under future climate change (Morgner et al., 2010;Nobrega and Grogan, 2007;Rogers et al., 2011;Schimel et al., 2004;Wahren et al., 2005;Welker et al., 2000).Apart from the spatial extent and duration of snowpack, the spatial heterogeneity of snow depth is an important factor in various terrestrial processes (Clark et al., 2011;Lundquist and Dettinger, 2005).The spatial distribution of snow not only affects the quantity of snowmelt discharge (Hartman et al., 1999;Luce et al., 1998), but also the water chemistry (Rohrbough et al., 2003;Wadham et al., 2006;Williams et al., 2001).Lawrence and Swenson (2011) demonstrated the importance of snow depth heterogeneity in predicting responses of the Arctic ecosystem to future climate change by performing idealized numerical simulations of shrub expansion across the pan-Arctic region using the Community Land Model (CLM4).Their results showed that an increase in active layer thickness (ALT) under shrubs was negated when spatial heterogeneity in snow cover due to wind driven snow redistribution was accounted for, resulting in an unchanged grid cell mean active layer thickness.López-Moreno et al. (2014) identified processes responsible for snow depth heterogeneity at three distinct spatial scales: microtopography at 1-10 m (Lopez-Moreno et al., 2011); wind induced lateral transport processes at 100-1000 m (Liston et al., 2007); and precipitation variability at catchment scales of 10 -1000 km (Sexstone and Fassnacht, 2014).
Large portions of the Arctic are characterized by polygonal ground features, which are formed in permafrost soil when frozen ground cracks due to thermal contraction during winter and ice wedges form within the upper several meters (Hinkel et al., 2005).
Polygons can be classified as 'low-centered' or 'high-centered' based on the relationship between their central and mean elevations.Polygonal ground features are dynamic components of the Arctic landscape in which the upper part of ice-wedge thaw under lowcentered polygon troughs leads to subsidence, eventually (~o(centuries)) converting the low-centered polygon into a high-centered polygon (Seppala et al., 1991).Microtopography of polygonal ground influences soil hydrologic and thermal conditions (Engstrom et al., 2005).In addition to controlling CO2 and CH4 emissions, soil moisture affects (1) partitioning of incoming radiation into latent, sensible, and ground heat fluxes (Hinzman and Kane, 1992;McFadden et al., 1998); (2) photosynthesis rates (McGuire et al., 2000;Oberbauer et al., 1991;Oechel et al., 1993;Zona et al., 2011); and (3) vegetation distributions (Wiggins, 1951).Our goals in this study include (1) analyzing the effects of spatially heterogeneous snow in polygonal ground on soil temperature and moisture and surface processes (e.g., surface energy budgets); (2) analyzing how model predictions are affected by inclusion of lateral subsurface hydrologic and thermal processes; and (3) developing and testing a three-dimensional version of the land model ALM (Tang and Riley, 2016;Zhu and Riley, 2015) integrated in the ACME Earth System Model (ESM).We note that the original version of ALM is equivalent to CLM4.5 (Koven et al., 2013;Oleson, 2013a), and represents vertical energy and water dynamics, including phase change.We expanded on that model to explicitly represent soil lateral energy and hydrological exchanges and fine-resolution snow redistribution (ALMv0-3D).We then applied ALMv0-3D to a transect across a polygonal tundra landscape at the Barrow Environmental Observatory in Alaska.After defining our study site, the model improvements, model tests against observations, and analyses, we apply the model to examine the effects of snow redistribution and lateral subsurface processes on snow micro-topographical heterogeneity, soil temperature, and the surface energy budget.

Study Area
Our analysis focuses on sites located near Barrow, Alaska (71.3 0 N, 156.5 0 W) from the long term Department of Energy (DOE) Next-Generation Ecosystem Experiment (NGEE-Arctic) project.The four primary NGEE-Arctic study sites (A, B, C, D) are located within the Barrow Environmental Observatory (BEO), which is situated on the Alaskan Coastal Plain.
The annual mean air temperature for our study sites is approximately -13°C (Walker et al., 2005) and mean annual precipitation is 106 mm with the majority of precipitation occurring during the summer season (Wu et al., 2013).The study site is underlain with continuous permafrost (Brown et al., 1980) and the annual maximum thaw depth (active layer depth) ranges between 30-90 cm (Hinkel et al., 2003).Although the overall topographic relief for the BEO is low, the four NGEE study sites have distinct

ALMv0 Description
We developed the capability to represent three-dimensional hydrology and thermal dynamics in ALMv0 (Zhu et al., 2016b), and call the new model ALMv0-3D.ALMv0 was derived from CLM4.5 (Ghimire et al., 2016;Koven et al., 2013), and is the land model integrated in the ACME Earth System Model (ESM).The model represents coupled plant biophysics, soil hydrology, and soil biogeochemistry (Oleson et al. 2013).We run ALMv0-3D here with prescribed plant phenology (called Satellite Phenology (SP) mode), since our focus is on the thermal dynamics of the system, rather than the C cycle dynamics.

Subsurface hydrology
The flow water in the unsaturated zone is given by the -based Richards equations as where [m of water m -3 of soil s -1 ] is volumetric sink of water.Darcy flux is given by where  [ms -1 ] is the hydraulic conductivity and  [m] is the soil matric potential.The hydraulic conductivity and soil matric potential are non-linear functions of volumetric soil moisture.ALMv0 uses the modified form of Richards equation of Zeng and Decker (2009) that computes Darcy flux as where C is a constant hydraulic potential above the water where  ![m] is the equilibrium soil matric potential.Substituting equations ( 3) and (4) into equation (1) yields the equation for the vertical transport of water in ALMv0: A finite volume spatial discretization and implicit temporal discretization with Taylor series expansion leads to a tri-diagonal system of equations.We extended this 1-D Richards equation to a 3-D representation integrated in ALMv0-3D, which is presented next.
We use a cell-centered finite volume discretization to decompose the spatial domain into  non-overlapping control volumes, Ω !, such that Ω = Ω !!!! ! and Γ !represents the boundary of the -th control volume.Applying a finite volume integral to equation (1) and the divergence theorem yields The spatially discretized equation for the -th grid cell that has  !volume and  !neighbors is given by For the sake of simplicity in presenting the discretized equation, we assume the 3-D grid is a Cartesian grid with each grid cell having a thickness of ∆, ∆, and ∆ in the , , and  directions, respectively.Using an implicit time integral, the 3-D discretized equation at time + 1 for a (, , ) control volume is given as where  !,  ! and  ! are Darcy flux in the , , and  directions, respectively and ∆ !,!,! !!! is the change in volumetric soil liquid water in time ∆.Using the same approach as Oleson The linearized Darcy fluxes in the  and  directions are computed similarly.Substituting 195 equation ( 9) in equation ( 8) results in a banded matrix of the form 196 where , , and  are subdiagonal entries; , , and  are superdiagonal entries;  is 197 diagonal entry of the banded matrix; and  is a column vector given by 198 The coefficients of equation ( 10) described in equation ( 11)-( 18) are for an internal grid cell with six neighbors.The coefficients for the top and bottom grid cells are modified for infiltration and interaction with the unconfined aquifer in the same manner as Oleson (2013b).Similarly, the coefficients for the grid cells on the lateral boundary are modified for a no-flux boundary condition.See Oleson (2013b) for details about the computation of hydraulic properties and derivative of Darcy flux with respect to soil liquid water content.

Subsurface thermal
ALMv0 solves a tightly coupled system of equations for soil, snow, and standing water temperature (Oleson, 2013a).The model solves the transient conservation of energy: where  is the volumetric heat capacity [J m -3 K -1 ], F is the heat flux [W m -2 ], and t is time The heat conduction flux is given by where  is thermal conductivity [W m -1 K -1 ] and T is temperature [K].Applying a finite volume integral to equation ( 20) and divergence theorem yields The spatially discretized equation for a -th grid cell that has  !volume and  !neighbors is given by Similar to the approach taken in Section 2.3.1,ALMv0-3D assumes a 3-D Cartesian grid with each grid cell having a thickness of ∆, ∆, and ∆ in the , , and  directions, respectively.Temporal integration of equation ( 22) is carried out using the Crank-Nicholson method that uses a linear combination of fluxes evaluated at time t and t + 1: where  is the weight in the Crank-Nicholson method and set to 0.5 in this study.218 Substituting a discretized form of heat flux using equation ( 20) in equation ( 23), results in 219 a banded matrix of the form 220 where , , and  are subdiagonal entries; , , and  are superdiagonal entries;  is 221 diagonal entry of the banded matrix; and  is a column vector given by 222 223 224 226 The coefficients of equation ( 24) described in equation ( 25)-( 32) are for an internal grid cell with six neighbors.The coefficients for the top and bottom grid cells are modified for presence of snow and/or standing water, and no-flux boundary.The coefficients for the grid cells on the lateral boundary are modified for a no-flux boundary condition.ALM handles ice-liquid phase transitions by first predicting temperatures at the end of a time step and then updating temperatures after accounting for deficits or excesses of energy during melting or freezing.See Oleson (2013b) for details about the computation of thermal properties and phase transition.

Numerical solution via PETSc
ALMv0, which considers flow only in the vertical direction, solves a tridiagonal and banded tridiagonal system of equations for water and energy transport, respectively.In ALMv0-3D, accounting for lateral flow in the subsurface results in a sparse linear system, equations ( 10) and ( 24), where the sparcity pattern of the linear system depends on grid cell connectivity.In this work, we use the PETSc (Portable, Extensible Toolkit for Scientific Computing) library (Balay et al., 2016) developed at the Argonne National Laboratory to solve the sparse linear systems.PETSc provides object-oriented data structures and solvers for scalable scientific computation on parallel supercomputers.

Snow Model and Redistribution
The snow model in ALMv0-3D is the same as that in the default ALMv0 and CLM4.5 (Anderson, 1976;Dai and Zeng, 1997;Jordan, 1991).The snow model allows for a dynamic snow depth and up to 5 snow layers, and explicitly solves the vertically-resolved mass and energy budgets.Snow aging, compaction, and phase change are all represented in the snow model formulation.Additionally, the snow model accounts for the influence of aerosols (including black and organic carbon and mineral dust) on snow radiative transfer (Oleson, 2013a).ALMv0 uses the methodology of Swenson and Lawrence (2012) to compute fractional snow cover area, which is appropriate for ESM-scale grid cells (~100 [km] x 100 [km]).Since the grid cell resolution in this work is sub-meter, we modified the fractional cover to be either 1 (when snow was present) or 0 (when snow was absent).Two main drivers of snow redistribution (SR) include topography and surface wind (Warscher et al., 2013); previous SR models include mechanistically- (Bartelt and Lehning, 2002;Liston and Elder, 2006) and empirically- (Frey and Holzmann, 2015;Helfricht et al., 2012) based approaches.To mimic the effects of wind, we used a conceptual model to simulated SR over the fine-resolution topography of our site by instantaneously re-distributing the incoming snow flux such that lower elevation areas (polygon center) receive snow before higher elevation areas (polygon rims).This relatively simple and parsimonious approach is reasonable given the observed snow depth heterogeneity, as described below, and small spatial extent of our domain.

System Characterization
Hydrologic and thermal properties differ by depth and landscape type.We used the horizontal distribution of OM organic matter from Wainwright et al. (2015) to infer soil hydrologic and thermal properties following the default representations in ALM.
Vegetation cover was classified as arctic shrubs in polygon centers and arctic grasses in polygon rims.The default representation of the plant wilting factor assigns a value of zero

Simulation Setup, Climate Forcing, and Analyses
Because of computational constraints, we investigated the role of snow redistribution and physics representation using a two-dimensional transect through site A (Figure 1 was used in the simulations.All simulations were performed in "SP" mode, i.e., Leaf Area Index (LAI) was prescribed from MODIS observations.
Simulations were run for 10 years using long-term climate data gathered at the Barrow, Alaska Observatory site (https://www.esrl.noaa.gov/gmd/obop/brw/)managed by the Global Monitoring Division of NOAA's Earth System Research Laboratory (Mefford et al., 1996).The missing precipitation time series was gap-filled using daily precipitation at the Barrow Regional Airport available from the Global Historical Climatology Network 3 Results and Discussion

Snow depth
In the absence of SR, predicted snow depth exactly follows the topography.With SR, a much larger dependence of winter-average snow depth on topography is predicted (Figure 2).Further, for the winter average, there are very small differences in snow depth between simulations with SR and 1D or 2D subsurface physics representations.Compared to observations, considering snow redistribution led to: (1) a factor of ~2 improvement in snow depth bias for the polygon center; (2) modest increase and decrease in average bias on the rims for September through February and March through June, respectively; and (3) a dramatic improvement in bias of the difference in snow depth between the polygon centers and rims (Figure 3).There was no discernible difference in snow depth bias between the 1D and 2D physics (Table 1), although the predicted subsurface temperature fields were different, as shown below.
The temporal variation of the mean snow depth (Figure 4a) and its spatial standard deviation (Figure 4b) also differed based on whether SR was considered, but was not affected by considering 2D thermal or hydrologic physics.With SR, the snow depth coefficient of variation (Figure 4c) was about 0.5 from December through the beginning of the snowmelt period, indicating relatively large spatial heterogeneity.Snapshots of simulated snow depth for the three simulation scenarios are included in Supplementary material (Error!Reference source not found.).

Soil Temperature and Active Layer Depth
Broadly, ALMv0-3D accurately predicted the polygon center soil temperature at depth intervals corresponding to the temperature probes (0-20 cm, 20-50 cm, 50-75 cm, and 75-100 cm; Figure 5a).Recall that the observed temperatures for the polygon center and rims were taken at single points in site A (Figure 1) while the predicted temperatures were calculated as averages across the transect for each of the two landscape position types.The model was able to simulate early freeze up of the soil column under the rims as compared to centers in November 2012 because of differences in accumulated snow pack.
The transition to thawed soil in the 0-20 cm depth interval in early June 2013 and the subsequent temperature dynamics over the summer were very well captured by ALMv0-3D.Minimum temperatures during the winter were also accurately predicted, although the temperatures in the deepest layer (75-100 cm) were overestimated by ~3°C in March.For figure clarity we did not indicate the standard deviation of the observations, but provide that information in Supplemental Material (Error!Reference source not found.-Error!Reference source not found.).
Similarly, the soil temperatures were accurately predicted in the polygon rims (Figure 5b).The largest discrepancies between measured and predicted soil temperatures were in the shallowest layer (0 -25 cm), where the predictions were up to a few °C cooler than some of the observations between December 2012 and March 2013.In the polygon center, a thicker snow pack acts as a heat insulator and keeps soil temperature higher in winter as compared to the polygon rims.Snow redistribution impacts spatial variability of soil temperature throughout the soil column.Absence of SR results in no significant spatial variability of soil temperature (Figure 6a).Inclusion of SR on the surface modifies the amount of energy exchanged between the snow and the top soil layer, thereby creating spatial variability in the temperature of the top soil, which propagates down into the soil column (Figure 6b).With SR, energy dissipation in the lateral direction reduces the penetration depth of the soil temperature spatial variance (compare Figure 6c and Figure 6b).
With 1D physics, the average spatial and temporal difference of the active layer depth (ALD) between simulations with and without SR was 1.7 cm (Figure 7a), and the absolute difference was 6.5 cm.As described above, we diagnosed the ALD to be the maximum soil depth during the summer at which vertically interpolated soil temperature is 0 °C.On average, the rims had ~10 cm shallower ALD with (blue line) than without (green line) SR, consistent with the loss of insulation from SR on the rims during the winter.In the centers (e.g., at location 42 -55 m), the thaw depth was deeper by ~5 cm with SR because of the higher snow depth there from SR.The effect of SR on the ALD was largest on the rims because, compared to centers, they (1) on average lost more snow with SR and (2) are more thermally conductive.Since rims are therefore colder at the time of snowmelt with SR, the ground heat flux during the subsequent summer was unable to thaw the soil column as deeply as when SR is ignored.For comparison, Atchley et al. (2015) found in their sensitivity analysis using the 1D version of ATS that SR resulted in deeper thaw depths in both polygon centers (by ~3 cm) and rims (~0.3 cm).Thus, there results for polygon centers are consistent in sign but lower in magnitude than ours, but opposite in sign for the rims.
Across ten years of simulation, the inter-annual variability (IAV) in ALD varied substantially between the three scenarios (Figure 7b).As expected, for the 1D physics without SR scenario (green line), the IAV in ALD was determined by landscape position because of differences in soil and vegetation parameters.With SR and 1D physics, the model shows largest differences over the rims, again highlighting the relatively larger effects of SR on the rim soil temperatures.The effect of 1D versus 2D physics on the ALD across the transect was modest (mean absolute difference ~3 cm).Generally, because 2D physics allows for lateral energy diffusion, the horizontal variation of ALD was slightly lower (i.e., the red line is smoother than the blue line; Figure 7a) than with 1D physics.This difference was also reflected in the thaw depth IAV across the transect, where 2D physics led to a smoother lateral profile of inter-annual variability than with 1D physics.
The impact of physics formulation (i.e., 1D or 2D) alone was investigated by analyzing differences between soil temperature profiles over time for polygon rims and centers in simulations with snow redistribution.Inclusion of 2D subsurface physics resulted in soil temperatures with depth and time that were lower in the polygon rims (Figure 8a) and higher in polygon centers (Figure 8b).Using the simulations from the scenario with SR and 2D physics, we evaluated the extent to which the soils under rims and centers can be separately considered as relatively homogeneous single column systems by evaluating the soil temperature standard deviation as a function of depth and time (Figure 9).During winter, both polygon rims and centers showed soil temperature spatial variability >1 °C up to a depth of ~2 [m].The soil temperature spatial variability in winter due to snow redistribution is dissipated over the summer.During the summer, polygon centers were relatively more homogeneous vertically compared to polygon rims.

Surface Energy Budget
Predicted monthly-and spatial-mean (μ) surface latent heat fluxes across the transect were very similar between the three scenarios (Figure 10a), with a growing  The predicted temporal monthly-mean and spatial-mean surface sensible heat fluxes across the transect were also similar between the three scenarios (Figure 11a), with a growing season mean absolute difference of < 3.5 W m -2 .Also, the sensible heat flux spatial variability differences occurred earlier than snowmelt, in contrast to the latent heat flux.Both the standard deviation and CV of the sensible heat fluxes were larger than those of the latent heat fluxes, with early season standard deviations of ~50 W m -2 (Figure 11b) and CV's of ~1.5 (Figure 11c).As for the latent heat fluxes, the differences in standard deviation and CV of sensible heat fluxes were small between the 1D and 2D scenarios with SR, arguing that the subsurface lateral energy exchanges associated with the 2D physics did not propagate to the mean surface heat fluxes.However, as for the latent heat flux, there was a relatively large difference in spatial variation between the scenarios with and without SR (e.g., of about 25 W m -2 in May; Figure 10b).

Soil Moisture
Neither SR nor 2D lateral physics affected the spatial mean moisture across time (not shown).However, the spatial heterogeneity of predicted soil moisture content differed substantially between scenarios during the snow free period (Figure 12).For the 1D simulations, the effect of SR was to increase the growing season soil moisture spatial heterogeneity by factors of 5.2 and 1.6 for 0-10 cm and 10-65 cm depth intervals, respectively (compare Figure 12a and Figure 12b).Compared to the 1D physics, simulating 2D thermal and hydrologic physics led to an overall reduction in the soil moisture spatial heterogeneity by factors of 0.8 and 0.7 for 0-10 cm and 10-65 cm depth intervals, respectively (compare Figure 12b and Figure 12c).Thus, with respect to dynamic spatial mean soil moisture, SR effects dominated those associated with lateral subsurface water movement .

Caveats and Future Work
The good agreement between ALMv0-3D predictions and soil temperature observations demonstrate the model's capabilities to represent this very spatially heterogeneous and complex system.However, several caveats to our conclusions remain due to uncertainties in model parameterizations, model structure, and climate forcing data.Because of computational constraints, we applied a 2D transect domain to the site, instead of a full 3D domain.We are working to improve the computational efficiency of the model, which will facilitate a thorough analysis of the effects of 3D subsurface energy and water fluxes.A related issue is our simplified treatment of surface water flows.A thorough analysis of the effects of surface water redistribution would require integration of a 2D surface thermal flow model with the ALMv0-3D in a 3D domain, which is another goal for our future work.However, we note that the good agreement using the 2D model domain supports the idea that a two-dimensional simplification may be appropriate for this system.
The expected geomorphological changes in these systems over the coming decades (e.g., Liljedahl et al. 2016), which will certainly affect soil temperature and moisture, are not currently represented in ALM, although incorporation of these processes is a long-term development goal.
The current representation of vegetation in ALMv0-3D for these polygonal tundra systems is over-simplified.For example, non-vascular plants (mosses and lichens) are not explicitly represented in the model, but can be responsible for a majority of evaporative losses (Miller et al., 1976) and are strongly influenced by near surface hydrologic conditions (Williams and Flanagan, 1996).Our use of the 'satellite phenology' mode, which imposes transient LAI profiles for each plant functional type in the domain, ignores the likely influence of nutrient constraints (Zhu et al., 2016a) on photosynthesis and therefore the surface energy budget.Other model simplifications, e.g., the simplified treatment of radiation competition may also be important, especially as simulations are extended over periods where vegetation change may occur (e.g., Grant 2016).

Summary and Conclusions
We analyzed the effects of microtopographical surface heterogeneity and lateral subsurface transport in a polygonal tundra landscape on soil temperature, soil moisture, and surface energy exchanges.Starting from the climate-scale land model ALMv0, we incorporated in ALMv0-3D numerical representations of subsurface water and energy lateral transport that are solved using PETSc.A simple method for redistributing incoming snow along the microtopographic transect was also integrated in the model.Overall, our analysis demonstrates the potential and value of explicitly representing 507 snow redistribution and lateral subsurface hydrologic and thermal dynamics in polygonal 508 ground systems and quantifies the effects of these processes on the resulting system states 509 and surface energy exchanges with the atmosphere.The integration of 3D subsurface 510 processes in the ACME Land Model also allows for a wide range of analyses heretofore 511  Table 1.Bias, root mean square error (RMSE), and correlation (R 2 ) between modeled and 515 observed snow depth at polygon center, rim and difference between center and rim for 516
microtopographic features: low-centered (A), high-centered (B), and transitional polygons (C, D).Contrasting polygon types are indicative of different stages of permafrost Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2017-71Manuscript under review for journal Geosci.Model Dev. Discussion started: 6 June 2017 c Author(s) 2017.CC BY 3.0 License.degradation and were the primary motivation behind the choice of study sites for the NGEE-Arctic project.LIDAR Digital Elevation Model (DEM) data were available at 0.25 m resolution for the region encompassing all four NGEE sites.In this work, we perform simulations along a two-dimensional transect in low-centered polygon Site-A as shown by the dotted line in Figure 1.
Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2017-71Manuscript under review for journal Geosci.Model Dev. Discussion started: 6 June 2017 c Author(s) 2017.CC BY 3.0 License.for a given soil layer when it's temperature falls below a threshold (Tthreshold) of -2 0 C.This default value leads to overly large predicted latent and sensible heat fluxes during winter, compared to nearby eddy covariance measurements.We modified Tthreshold to be 0 0 C in this study, resulting in improved predicted wintertime latent heat fluxes compared to the default version of the model (Error!Reference source not found.).Although biases compared to the observations remain, particularly for sensible heat fluxes in the spring, the improvement is substantial and, given the observational uncertainties, we believe sufficient to justify our use of the model for investigations of the role of snow heterogeneity in this polygonal tundra system.
).The transect was 104[m]  long and 45[m]  deep that was discretized horizontally with a grid spacing of 0.25 [m] and an exponentially varying layer thickness in the vertical with 30 soil layers.No flow conditions for mass and energy were imposed on the east, west, and bottom boundaries of the domain.Temporal discretization of 30[min]

(
http://www1.ncdc.noaa.gov/pub/data/ghcn/daily). We tested the model by comparing predictions to high-frequency observations of snow depth and vertically resolved soil temperature for September 2012 -September 2013.Temperature observations were taken at discrete locations in a polygon center and rim (Figure1), and were combined to analyze comparable landscape positions in the simulations (Figure2).After testing, the model was used to investigate the effect of snow redistribution and 2D subsurface hydrologic and thermal physics by analyzing three scenarios: (1) no snow Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2017-71Manuscript under review for journal Geosci.Model Dev. Discussion started: 6 June 2017 c Author(s) 2017.CC BY 3.0 License.redistribution and 1D physics; (2) snow redistribution and 1D physics; and (3) snow redistribution and 2D physics.Between these scenarios, we compared vertically-resolved soil temperature and liquid saturation, active layer depth, and mean and spatial variation of latent and sensible heat fluxes across the 10 years of simulations.For each soil column, the simulated soil temperature was interpolated vertically and the active layer depth was estimated as the maximum depth that had above-freezing soil temperature.
Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2017-71Manuscript under review for journal Geosci.Model Dev. Discussion started: 6 June 2017 c Author(s) 2017.CC BY 3.0 License. in this work or that of Kumar et al. (2016), while the ATS parameterizations were tuned to match the soil temperature profile.
season mean difference of < 1.0 [W m-2 ].However, the spatial variability (SV = σ; Figure10b) and coefficient of variation (CV = σ/μ; Figure10c) of latent heat fluxes were different between the scenarios with SR (1D and 2D physics) and without SR.With SR, the latent heat flux spatial standard deviation peaked after snowmelt and declined until the fall when snow began, from about ~100% to 10% of the mean.This relatively larger spatial variation in latent heat flux occurred because of large spatial heterogeneity in near surface soil moisture in the beginning of summer, indicating a residual effect of SR from the previous winter.
Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2017-71Manuscript under review for journal Geosci.Model Dev. Discussion started: 6 June 2017 c Author(s) 2017.CC BY 3.0 License.Over the observational record, ALMv0-3D with snow redistribution and lateral heat and hydrological fluxes accurately predicted snow depth and soil temperature vertical profiles in the polygon rims and centers (overall bias, RMSE, and R 2 of 0.59 0 C, 1.82 0 C and 0.99, respectively).In the rims, the transition to thawed soil in spring, summer temperature dynamics, and minimum temperatures during the winter were all accurately predicted.In the centers, a ~2°C warm bias in April in the 75-100 cm soil layer was predicted, although this bias disappeared during snowmelt.The spatial heterogeneity of snow depth during the winter due to snow redistribution generated surface soil temperature heterogeneity that propagated into the soil over time.The temporal and spatial variation of snow depth was affected by snow redistribution, but not by lateral thermal and hydrologic transport.Both snow redistribution and lateral thermal fluxes affected spatial variability of soil temperatures.Energy dissipation in the lateral direction reduced the depth to which soil temperature variance penetrated.Snow redistribution led to ~10 cm shallower active layer depths under the polygon rims because of the residual effect of reduced insulation during the winter.In contrast, snow redistribution led to ~5 cm deeper active layers under the polygon centers.The effect of lateral energy fluxes on active layer depths was ~3 cm.Compared to 1D physics, the 2D subsurface physics led to lower (higher) soil temperatures with depth and time in the polygon rims (centers).The larger than 1 °C wintertime spatial temperature variability down to ~2 m depth in rims and centers indicates the uncertainty associated with considering rims and centers as separate 1D columns.During the summer, polygon center temperatures were relatively more vertically homogeneous than temperatures in the rims.The monthly-and spatial-mean predicted latent and sensible heat fluxes were unaffected by snow redistribution and lateral heat and hydrological fluxes.However, snow redistribution led to spatial heterogeneity in surface energy fluxes and soil moisture during the summer.Excluding lateral subsurface hydrologic and thermal processes led to an over prediction of spatial variability in soil moisture and soil temperature because subsurface gradients were artificially prevented from laterally dissipating over time.Snow redistribution effects on soil moisture heterogeneity were larger than those associated with lateral thermal fluxes.

Figure 1 Figure 2 .Figure 3 Figure
Figure 1 The NGEE-Arctic study area A, which characterized as a low-centered polygon field.Dotted line indicate the transect along which simulation in this paper are preformed to demonstrate the effects of snow redistribution on soil temperature.The locations where snow and temperature sensors are installed within the study site are denoted by triangle and circle, respectively.

Figure 5
Figure 5 Comparison of soil temperature observations and predictions in polygon centers 554 (a) and rims (b).Simulation was performed with snow redistribution on and 2D subsurface 555 physics, between September 2012 and September 2013.Simulation results are shown at an 556 interval of 10 days, while observations are shown at daily interval 557

Figure 6 569 Figure 7 Figure 8 Figure 9
Figure 6 Simulated daily spatial standard deviation averaged across 10-year of near surface soil temperature for simulation performed with snow redistribution turned off and 1D subsurface physics (top panel); snow redistribution turned on and 1D subsurface physics (middle panel); and snow redistribution turned on and 2D subsurface physics (bottom panel).