Impacts of microtopographic snow redistribution and lateral subsurface processes on hydrologic and thermal states in an Arctic polygonal ground ecosystem: a case study using ELM-3D v1.0

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 E3SM to redistribute incoming snow by accounting for microtopography and incorporated subsurface lateral transport of water and energy (ELM-3D v1.0). Multiple 10-year-long simulations were performed for a transect across a 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 predictions better agreed (higher R², lower bias and RMSE) with 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² 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 maximum thaw depths was modest, with mean absolute differences of ∼ 3 cm. Our integration of three-dimensional subsurface hydrologic and thermal subsurface dynamics in the E3SM land model will facilitate a wide range of analyses heretofore impossible in an ESM context.

1 Introduction

The northern circumpolar permafrost region, which contains ∼  1700 Pg of organic carbon down to 3 m (Tarnocai et al., 2009), is predicted to experience disproportionately larger future warming compared to the tropics and temperate latitudes (Holland and Bitz, 2003). Recent warming in the Arctic has led to changes in lake area (Smith et al., 2005), snow cover duration and extent (Callaghan et al., 2011a), vegetation cover (Sturm et al., 2005), growing season length (Smith et al., 2004), thaw depth (Schuur et al., 2008), permafrost stability (Jorgenson et al., 2006), and land–atmosphere feedbacks (Euskirchen et al., 2009). Future predictions of Arctic warming include northward expansion of shrub cover in tundra (Sturm et al., 2001; Tape et al., 2006), decreases in snow cover duration (Callaghan et al., 2011a), and emissions of CO₂ and CH₄ from decomposition of belowground soil organic matter (Koven et al., 2011; Schaefer et al., 2011; Schuur and Abbott, 2011; Xu et al., 2016).

Several recent modeling studies have predicted a positive global 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., 2005, 2007; 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., 2007; Randerson et al., 2015).

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 meltwater). 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 CO₂ 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). As synthesized by López-Moreno et al. (2015), the following processes are responsible for snow depth heterogeneity at three distinct spatial scales: microtopography at 1–10 m (López-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). 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, which is the maximum annual thaw depth, 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.

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 low-centered 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 CO₂ and CH₄ 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 the inclusion of lateral subsurface hydrologic and thermal processes; and (3) developing and testing a three-dimensional version of the E3SM Land Model (ELM; Tang and Riley, 2016; Zhu and Riley, 2015), called ELM-3D v1.0 (hereafter ELM-3D). We then applied ELM-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 microtopographical heterogeneity, soil temperature, and the surface energy budget.

2 Methodology

2.1 Study area

Our analysis focuses on sites located near Barrow, Alaska (71.3^∘ N, 156.5^∘ 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 (Sellmann et al., 1975) and the annual maximum thaw depth (active layer depth) ranges between 30 and 90 cm (Hinkel et al., 2003). Although the overall topographic relief for the BEO is low, the four NGEE study sites have distinct microtopographic features: low-centered (A), high-centered (B), and transitional polygons (C, D). Contrasting polygon types are indicative of different stages of permafrost 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 Fig. 1.

Figure 1 The NGEE-Arctic study area A, which was characterized as a low-centered polygon field. The dotted line indicates the transect along which the simulations in this paper are performed 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 triangles and circles, respectively.

2.2 ELMv0 description

The original version of ELM is equivalent to CLM4.5 (Ghimire et al., 2016; Koven et al., 2013; Oleson et al., 2013) and represents vertical energy and water dynamics, including phase change. We developed ELM-3D by expanding on that model to explicitly represent soil lateral energy and hydrological exchanges and fine-resolution snow redistribution. We run ELM-3D here with prescribed plant phenology (called the “satellite phenology” (SP) mode) since our focus is on thermal dynamics of the system, rather than C cycle dynamics.

2.3 Representing two-dimensional and three-dimensional physics

2.3.1 Subsurface hydrology

The flow of water in the unsaturated zone is given by the θ-based Richards equation as

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial \mathit{\theta }}{\partial t}}=-\mathrm{\nabla }\cdot \mathit{q}-Q,\end{array}$$

where θ (m³ m⁻³) is the volumetric soil water content, t (s) is time, q (ms⁻¹) is the Darcy flux, and Q (m⁻³ of water m⁻³ of soil s⁻¹) is the volumetric sink of water. The Darcy flux is given by

$$\begin{array}{}\text{(2)}& \mathit{q}=-k\mathrm{\nabla }(\mathit{\psi }+z),\end{array}$$

where k (ms⁻¹) is the hydraulic conductivity, ψ (m) is the soil matric potential, and z (m) is the height above a reference datum. The hydraulic conductivity and soil matric potential are nonlinear functions of volumetric soil moisture. ELMv0 uses the modified form of the Richards equation of Zeng and Decker (2009) that computes the Darcy flux as

$$\begin{array}{}\text{(3)}& \mathit{q}=-k\mathrm{\nabla }(\mathit{\psi }+z-C),\end{array}$$

where C is a constant hydraulic potential above the water table, z_∇, given as

$$\begin{array}{}\text{(4)}& C={\mathit{\psi }}_{\mathrm{E}}+z={\mathit{\psi }}_{\mathrm{sat}}{\left[{\displaystyle \frac{{\mathit{\theta }}_{\mathrm{E}}\left(z\right)}{{\mathit{\theta }}_{\mathrm{sat}}}}\right]}^{-B}+z={\mathit{\psi }}_{\mathrm{sat}}+z_{\mathrm{\nabla }},\end{array}$$

where ψ_E (m) is the equilibrium soil matric potential, ψ_sat (m) is the saturated soil matric potential, θ_E (m³ m⁻³) is the volumetric soil water content at equilibrium soil matric potential, θ_sat (m³ m⁻³) is the volumetric soil water content at saturation, z_∇ (m) is the height of water table above the reference datum, and B (−) is a fitting parameter for soil water characteristic curves. Substituting Eqs. (3) and (4) into Eq. (1) yields the equation for the vertical transport of water in ELMv0:

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\partial \mathit{\theta }}{\partial t}}={\displaystyle \frac{\partial }{\partial z}}\left[k\left({\displaystyle \frac{\partial \left(\mathit{\psi }-{\mathit{\psi }}_{\mathrm{E}}\right)}{\partial z}}\right)\right]-Q.\end{array}$$

A finite volume spatial discretization and implicit temporal discretization with Taylor series expansion leads to a tridiagonal system of equations. We extended this 1-D (one-dimensional) Richards equation to a 3-D (three-dimensional) representation integrated in ELM-3D, which is presented next.

We use a cell-centered finite volume discretization to decompose the spatial domain into N non-overlapping control volumes, Ω_n, such that $\mathrm{\Omega }={\bigcup }_{n=\mathrm{1}}^N{\mathrm{\Omega }}_n$, and Γ_n represents the boundary of the nth control volume. Applying a finite volume integral to Eq. (1) and the divergence theorem yields

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\partial }{\partial t}}\underset{{\mathrm{\Omega }}_n}{\int }\mathit{\theta }\mathrm{d}V=-\underset{{\mathrm{\Gamma }}_n}{\int }\left(\mathit{q}\cdot \mathrm{d}\mathit{A}\right)-\underset{{\mathrm{\Omega }}_{\mathrm{n}}}{\int }Q\mathrm{d}V.\end{array}$$

The spatially discretized equation for the nth grid cell that has V_n volume and n^′ neighbors is given by

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\mathrm{d}{\mathit{\theta }}_n}{\mathrm{d}t}}{V}_n=-\sum _{n^{\prime }}\left({\mathit{q}}_{n{n}^{\prime }}\cdot {\mathit{A}}_{n{n}^{\prime }}\right)-Q{V}_n.\end{array}$$

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 Δx, Δy, and Δz in the x, y, and z directions, respectively. Using an implicit time integral, the 3-D discretized equation at time t+1 for a ($i,j,k)$ control volume is given as

$$\begin{array}{ll}\text{(8)}& {\displaystyle }\left({\displaystyle \frac{\mathrm{\Delta }{\mathit{\theta }}_{i,j,k}^{t+\mathrm{1}}}{\mathrm{\Delta }t}}\right)V_{i,j,k}& {\displaystyle }=\left(q_{x_{i-\mathrm{1}/\mathrm{2},j,k}}^{t+\mathrm{1}}-q_{x_{i+\mathrm{1}/\mathrm{2},j,k}}^{t+\mathrm{1}}\right)\mathrm{\Delta }y\mathrm{\Delta }z{\displaystyle }& {\displaystyle }+\left(q_{y_{i,j-\mathrm{1}/\mathrm{2},k}}^{t+\mathrm{1}}-q_{y_{i,j+\mathrm{1}/\mathrm{2},k}}^{t+\mathrm{1}}\right)\mathrm{\Delta }x\mathrm{\Delta }z\\ {\displaystyle }& {\displaystyle }+\left(q_{z_{i,j,k-\mathrm{1}/\mathrm{2}}}^{t+\mathrm{1}}-q_{z_{i,j,k+\mathrm{1}/\mathrm{2}}}^{t+\mathrm{1}}\right)\mathrm{\Delta }x\mathrm{\Delta }y-Q{V}_{i,j,k},\end{array}$$

where q_x, q_y, and q_z are the Darcy flux in the x, y, and z directions, respectively, and $\mathrm{\Delta }{\mathit{\theta }}_{i,j,k}^{t+\mathrm{1}}$ is the change in volumetric soil liquid water in time, Δt. Using the same approach as Oleson et al. (2013), the Darcy flux in all three directions is linearized about θ using Taylor series expansion. The linearized Darcy flux in the x direction at the ($i-\mathrm{1}/\mathrm{2},j,k)$ interface is a function of ${\mathit{\theta }}_{i-\mathrm{1},j,k}$ and ${\mathit{\theta }}_{i,j,k}$:

$$\begin{array}{ll}\text{(9)}& {\displaystyle }{q}_{x_{i-\mathrm{1}/\mathrm{2},j,k}}^{t+\mathrm{1}}& {\displaystyle }=q_{x_{i-\mathrm{1}/\mathrm{2},j,k}}^t+{\displaystyle \frac{\partial q_{x_{i-\mathrm{1}/\mathrm{2},j,k}}^t}{\partial {\mathit{\theta }}_{i-\mathrm{1},j,k}}}\mathrm{\Delta }{\mathit{\theta }}_{i-\mathrm{1},j,k}^{t+\mathrm{1}}{\displaystyle }& {\displaystyle }+{\displaystyle \frac{\partial q_{x_{i-\mathrm{1}/\mathrm{2},j,k}}^t}{\partial {\mathit{\theta }}_{i,j,k}}}\mathrm{\Delta }{\mathit{\theta }}_{i+\mathrm{1},j,k}^{t+\mathrm{1}}.\end{array}$$

The linearized Darcy fluxes in the y and z directions are computed similarly. Substituting Eq. (9) in Eq. (8) results in a banded matrix of the form

$$\begin{array}{ll}\text{(10)}& {\displaystyle }\mathit{\alpha }\mathrm{\Delta }{\mathit{\theta }}_{i-\mathrm{1},j,k}^{t+\mathrm{1}}& {\displaystyle }+\mathit{\beta }\mathrm{\Delta }{\mathit{\theta }}_{i,j-\mathrm{1},k}^{t+\mathrm{1}}+\mathit{\gamma }\mathrm{\Delta }{\mathit{\theta }}_{i,j,k-\mathrm{1}}^{t+\mathrm{1}}+\mathit{\eta }\mathrm{\Delta }{\mathit{\theta }}_{i+\mathrm{1},j,k}^{t+\mathrm{1}}{\displaystyle }& {\displaystyle }+\mathit{\mu }\mathrm{\Delta }{\mathit{\theta }}_{i,j+\mathrm{1},k}^{t+\mathrm{1}}+\mathit{\varphi }\mathrm{\Delta }{\mathit{\theta }}_{i,j,k+\mathrm{1}}^{t+\mathrm{1}}+\mathit{\zeta }\mathrm{\Delta }{\mathit{\theta }}_{i,j,k}^{t+\mathrm{1}}=\mathit{\phi },\end{array}$$

where α, β, and γ are subdiagonal entries; η, μ, and ϕ are superdiagonal entries; ζ is a diagonal entry of the banded matrix and is given by

$$\begin{array}{ll}\text{(11)}& {\displaystyle }& {\displaystyle }\mathit{\alpha }={\displaystyle \frac{\partial q_{x_{i-\mathrm{1}/\mathrm{2},j,k}}^t}{\partial {\mathit{\theta }}_{i-\mathrm{1},j,k}}}\mathrm{\Delta }y\mathrm{\Delta }z\text{(12)}& {\displaystyle }& {\displaystyle }\mathit{\beta }={\displaystyle \frac{\partial q_{y_{i,j-\mathrm{1}/\mathrm{2},k}}^t}{\partial {\mathit{\theta }}_{i,j-\mathrm{1},k}}}\mathrm{\Delta }x\mathrm{\Delta }z\text{(13)}& {\displaystyle }& {\displaystyle }\mathit{\gamma }={\displaystyle \frac{\partial q_{z_{i,j,k-\mathrm{1}/\mathrm{2}}}^t}{\partial {\mathit{\theta }}_{i,j,k-\mathrm{1}}}}\mathrm{\Delta }x\mathrm{\Delta }y\text{(14)}& {\displaystyle }& {\displaystyle }\mathit{\eta }={\displaystyle \frac{\partial q_{x_{i+\mathrm{1}/\mathrm{2},j,k}}^t}{\partial {\mathit{\theta }}_{i+\mathrm{1},j,k}}}\mathrm{\Delta }y\mathrm{\Delta }z\text{(15)}& {\displaystyle }& {\displaystyle }\mathit{\mu }={\displaystyle \frac{\partial q_{y_{i,j+\mathrm{1}/\mathrm{2},k}}^t}{\partial {\mathit{\theta }}_{i,j+\mathrm{1},k}}}\mathrm{\Delta }x\mathrm{\Delta }z\text{(16)}& {\displaystyle }& {\displaystyle }\mathit{\varphi }={\displaystyle \frac{\partial q_{z_{i,j,k+\mathrm{1}/\mathrm{2}}}^t}{\partial {\mathit{\theta }}_{i,j,k+\mathrm{1}}}}\mathrm{\Delta }x\mathrm{\Delta }y\text{(17)}& {\displaystyle }& {\displaystyle }\mathit{\zeta }=\left({\displaystyle \frac{\partial q_{x_{i-\mathrm{1}/\mathrm{2},j,k}}^t}{\partial {\mathit{\theta }}_{i,j,k}}}-{\displaystyle \frac{\partial q_{x_{i+\mathrm{1}/\mathrm{2},j,k}^t}}{\partial {\mathit{\theta }}_{i,j,k}}}\right)\mathrm{\Delta }y\mathrm{\Delta }z{\displaystyle }& {\displaystyle }+\left({\displaystyle \frac{\partial q_{y_{i,j-\mathrm{1}/\mathrm{2},k}}^t}{\partial {\mathit{\theta }}_{i,j,k}}}-{\displaystyle \frac{\partial q_{y_{i,j+\mathrm{1}/\mathrm{2},k}^t}}{\partial {\mathit{\theta }}_{i,j,k}}}\right)\mathrm{\Delta }x\mathrm{\Delta }z\\ {\displaystyle }& {\displaystyle }+\left({\displaystyle \frac{\partial q_{z_{i,j,k-\mathrm{1}/\mathrm{2}}}^t}{\partial {\mathit{\theta }}_{i,j,k}}}-{\displaystyle \frac{\partial q_{z_{i,j,k+\mathrm{1}/\mathrm{2}}}^t}{\partial {\mathit{\theta }}_{i,j,k}}}\right)\mathrm{\Delta }x\mathrm{\Delta }y-{\displaystyle \frac{\mathrm{\Delta }x\mathrm{\Delta }x\mathrm{\Delta }z}{\mathrm{\Delta }t}}.\end{array}$$

The column vector φ is given by

$$\begin{array}{ll}\text{(18)}& {\displaystyle }\mathit{\phi }=& {\displaystyle }-\left(q_{x_{i-\frac{\mathrm{1}}{\mathrm{2}},j,k}}^t-q_{x_{i+\frac{\mathrm{1}}{\mathrm{2}},j,k}}^t\right)\mathrm{\Delta }y\mathrm{\Delta }z{\displaystyle }& {\displaystyle }-\left(q_{y_{i,j-\frac{\mathrm{1}}{\mathrm{2}},k}}^t-q_{y_{i,j+\frac{\mathrm{1}}{\mathrm{2}},k}^t}\right)\mathrm{\Delta }x\mathrm{\Delta }z\\ {\displaystyle }& {\displaystyle }-\left(q_{z_{i,j,k-\frac{\mathrm{1}}{\mathrm{2}}}}^t-q_{z_{i,j,k+\frac{\mathrm{1}}{\mathrm{2}}}}^t\right)\mathrm{\Delta }x\mathrm{\Delta }y+Q_{i,j,k}^{t+\mathrm{1}}\mathrm{\Delta }x\mathrm{\Delta }x\mathrm{\Delta }z.\end{array}$$

The coefficients of Eq. (10) described in Eqs. (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 et al. (2013). Similarly, the coefficients for the grid cells on the lateral boundary are modified for a no-flux boundary condition. See Oleson et al. (2013) for details about the computation of hydraulic properties and derivative of the Darcy flux with respect to soil liquid water content.

2.3.2 Subsurface thermal

ELMv0 solves a tightly coupled system of equations for soil, snow, and standing water temperature (Oleson et al., 2013). The model solves the transient conservation of energy:

$$\begin{array}{}\text{(19)}& c{\displaystyle \frac{\partial T}{\partial t}}=-\mathrm{\nabla }\cdot \mathbf{F},\end{array}$$

where c is the volumetric heat capacity (J m⁻³ K⁻¹), F is the heat flux (W m⁻²), and t is time (s). The heat conduction flux is given by

$$\begin{array}{}\text{(20)}& \mathbf{F}=-\mathit{\lambda }\mathrm{\nabla }T,\end{array}$$

where λ is thermal conductivity (W m⁻¹ K⁻¹) and T is temperature (K). Applying a finite volume integral to Eq. (20) and divergence theorem yields

$$\begin{array}{}\text{(21)}& c{\displaystyle \frac{\partial }{\partial t}}\underset{{\mathrm{\Omega }}_n}{\int }T=-\underset{{\mathrm{\Gamma }}_{\mathrm{n}}}{\int }\mathit{F}\cdot \mathrm{d}\mathit{A}.\end{array}$$

The spatially discretized equation for a nth grid cell that has V_n volume and n^′ neighbors is given by

$$\begin{array}{}\text{(22)}& c_n{\displaystyle \frac{\mathrm{d}{T}_n}{\mathrm{d}t}}{V}_n=-\sum _{n^{\prime }}\left({\mathit{F}}_{n{n}^{\prime }}\cdot {\mathit{A}}_{n{n}^{\prime }}\right).\end{array}$$

Similar to the approach taken in Sect. 2.3.1, ELM-3D assumes a 3-D Cartesian grid, with each grid cell having a thickness of Δx, Δy, and Δz in the x, y, and z directions, respectively. Temporal integration of Eq. (22) is carried out using the Crank–Nicholson method that uses a linear combination of fluxes evaluated at time t and t+1:

$$\begin{array}{ll}\text{(23)}& {\displaystyle }& {\displaystyle }{c}_{n_{i,j,k}}{\displaystyle \frac{\left(T_{i,j,k}^{t+\mathrm{1}}-T_{i,j,k}^t\right)}{\mathrm{\Delta }t}}\mathrm{\Delta }x\mathrm{\Delta }y\mathrm{\Delta }z{\displaystyle }& {\displaystyle }=\mathit{\omega }\left\{\left(F_{x_{i-\frac{\mathrm{1}}{\mathrm{2}},j,k}}^t-F_{x_{i+\frac{\mathrm{1}}{\mathrm{2}},j,k}}^t\right)\mathrm{\Delta }y\mathrm{\Delta }z\right.\\ {\displaystyle }& {\displaystyle }+\left(F_{y_{i,j-\frac{\mathrm{1}}{\mathrm{2}},k}}^t-F_{y_{i,j+\frac{\mathrm{1}}{\mathrm{2}},k}}^t\right)\mathrm{\Delta }x\mathrm{\Delta }z\\ {\displaystyle }& {\displaystyle }\left.+\left(F_{z_{i,j,k-\frac{\mathrm{1}}{\mathrm{2}}}}^t-F_{z_{i,j,k+\frac{\mathrm{1}}{\mathrm{2}}}}^t\right)\mathrm{\Delta }x\mathrm{\Delta }y\right\}\\ {\displaystyle }& {\displaystyle }+\left(\mathrm{1}-\mathit{\omega }\right)\left\{\left(F_{x_{i-\frac{\mathrm{1}}{\mathrm{2}},j,k}}^{t+\mathrm{1}}-F_{x_{i+\frac{\mathrm{1}}{\mathrm{2}},j,k}}^{t+\mathrm{1}}\right)\mathrm{\Delta }y\mathrm{\Delta }z\right.\\ {\displaystyle }& {\displaystyle }+\left(F_{y_{i,j-\frac{\mathrm{1}}{\mathrm{2}},k}}^{t+\mathrm{1}}-F_{y_{i,j+\frac{\mathrm{1}}{\mathrm{2}},k}}^{t+\mathrm{1}}\right)\mathrm{\Delta }x\mathrm{\Delta }z\\ {\displaystyle }& {\displaystyle }\left.+\left(F_{z_{i,j,k-\frac{\mathrm{1}}{\mathrm{2}}}}^{t+\mathrm{1}}-F_{z_{i,j,k+\frac{\mathrm{1}}{\mathrm{2}}}}^t+\mathrm{1}\right)\mathrm{\Delta }x\mathrm{\Delta }y\right\},\end{array}$$

where ω is the weight in the Crank–Nicholson method and is set to 0.5 in this study. Substituting a discretized form of heat flux using Eq. (20) in Eq. (23) results in a banded matrix of the form

$$\begin{array}{ll}\text{(24)}& {\displaystyle }\mathit{\alpha }{T}_{i-\mathrm{1},j,k}^{t+\mathrm{1}}& {\displaystyle }+\mathit{\beta }{T}_{i,j-\mathrm{1},k}^{t+\mathrm{1}}+\mathit{\gamma }{T}_{i,j,k-\mathrm{1}}^{t+\mathrm{1}}+\mathit{\eta }{T}_{i+\mathrm{1},j,k}^{t+\mathrm{1}}{\displaystyle }& {\displaystyle }+\mathit{\mu }{T}_{i,j+\mathrm{1},k}^{t+\mathrm{1}}++\mathit{\varphi }{T}_{i,j,k+\mathrm{1}}^{t+\mathrm{1}}+\mathit{\zeta }\mathrm{\Delta }{T}_{i,j,k}^{t+\mathrm{1}}=\mathit{\phi },\end{array}$$

where α, β, and γ are subdiagonal entries; η, μ, and ϕ are superdiagonal entries; ζ is a diagonal entry of the banded matrix and is given by

$$\begin{array}{ll}\text{(25)}& {\displaystyle }& {\displaystyle }\mathit{\alpha }=\left({\displaystyle \frac{-\left(\mathrm{1}-\mathit{\omega }\right)\mathrm{\Delta }t}{c_{n_{i,j,k}}\mathrm{\Delta }x}}\right)\left({\displaystyle \frac{{\mathit{\lambda }}_{i-\mathrm{1}/\mathrm{2},j,k}}{x_{i,j,k}-x_{i-\mathrm{1},j,k}}}\right)\text{(26)}& {\displaystyle }& {\displaystyle }\mathit{\beta }=\left({\displaystyle \frac{-\left(\mathrm{1}-\mathit{\omega }\right)\mathrm{\Delta }t}{c_{n_{i,j,k}}\mathrm{\Delta }y}}\right)\left({\displaystyle \frac{{\mathit{\lambda }}_{i,j-\mathrm{1}/\mathrm{2},k}}{y_{i,j,k}-y_{i-\mathrm{1},j,k}}}\right)\text{(27)}& {\displaystyle }& {\displaystyle }\mathit{\gamma }=\left({\displaystyle \frac{-\left(\mathrm{1}-\mathit{\omega }\right)\mathrm{\Delta }t}{c_{n_{i,j,k}}\mathrm{\Delta }z}}\right)\left({\displaystyle \frac{{\mathit{\lambda }}_{i,j,k-\mathrm{1}/\mathrm{2}}}{z_{i,j,k}-z_{i,j,k-\mathrm{1}}}}\right)\text{(28)}& {\displaystyle }& {\displaystyle }\mathit{\eta }=\left({\displaystyle \frac{-\left(\mathrm{1}-\mathit{\omega }\right)\mathrm{\Delta }t}{c_{n_{i,j,k}}\mathrm{\Delta }x}}\right)\left({\displaystyle \frac{{\mathit{\lambda }}_{i+\mathrm{1}/\mathrm{2},j,k}}{x_{i+\mathrm{1},j,k}-x_{i,j,k}}}\right)\text{(29)}& {\displaystyle }& {\displaystyle }\mathit{\mu }=\left({\displaystyle \frac{-\left(\mathrm{1}-\mathit{\omega }\right)\mathrm{\Delta }t}{c_{n_{i,j,k}}\mathrm{\Delta }y}}\right)\left({\displaystyle \frac{{\mathit{\lambda }}_{i-\mathrm{1}/\mathrm{2},j,k}}{y_{i+\mathrm{1},j,k}-y_{i,j,k}}}\right)\text{(30)}& {\displaystyle }& {\displaystyle }\mathit{\varphi }=\left({\displaystyle \frac{-\left(\mathrm{1}-\mathit{\omega }\right)\mathrm{\Delta }t}{c_{n_{i,j,k}}\mathrm{\Delta }z}}\right)\left({\displaystyle \frac{{\mathit{\lambda }}_{i-\mathrm{1}/\mathrm{2},j,k}}{z_{i+\mathrm{1},j,k}-z_{i,j,k}}}\right)\text{(31)}& {\displaystyle }& {\displaystyle }\mathit{\zeta }=\mathrm{1}+\left({\displaystyle \frac{\left(\mathrm{1}-\mathit{\omega }\right)\mathrm{\Delta }t}{c_{n_{i,j,k}}\mathrm{\Delta }x}}\right)\left[{\displaystyle \frac{{\mathit{\lambda }}_{i-\mathrm{1}/\mathrm{2},j,k}}{x_{i,j,k}-x_{i-\mathrm{1},j,k}}}\right.{\displaystyle }& {\displaystyle }\left.+{\displaystyle \frac{{\mathit{\lambda }}_{i+\mathrm{1}/\mathrm{2},j,k}}{x_{i+\mathrm{1},j,k}-x_{i,j,k}}}\right]+\left({\displaystyle \frac{\left(\mathrm{1}-\mathit{\omega }\right)\mathrm{\Delta }t}{c_{n_{i,j,k}}\mathrm{\Delta }y}}\right)\\ {\displaystyle }& {\displaystyle }\left[{\displaystyle \frac{{\mathit{\lambda }}_{i,j-\mathrm{1}/\mathrm{2},k}}{y_{i,j,k}-y_{i-\mathrm{1},j,k}}}+{\displaystyle \frac{{\mathit{\lambda }}_{i-\mathrm{1}/\mathrm{2},j,k}}{y_{i+\mathrm{1},j,k}-y_{i,j,k}}}\right]\\ {\displaystyle }& {\displaystyle }+\left({\displaystyle \frac{\left(\mathrm{1}-\mathit{\omega }\right)\mathrm{\Delta }t}{c_{n_{i,j,k}}\mathrm{\Delta }z}}\right)\left[{\displaystyle \frac{{\mathit{\lambda }}_{i,j,k-\mathrm{1}/\mathrm{2}}}{z_{i,j,k}-z_{i,j,k-\mathrm{1}}}}+{\displaystyle \frac{{\mathit{\lambda }}_{i-\mathrm{1}/\mathrm{2},j,k}}{z_{i+\mathrm{1},j,k}-z_{i,j,k}}}\right].\end{array}$$

The column vector φ is given by

$$\begin{array}{ll}\text{(32)}& {\displaystyle }\mathit{\phi }& {\displaystyle }=T_{i,j,k}^t+\left({\displaystyle \frac{\mathit{\omega }\mathrm{\Delta }t}{c_{n_{i,j,k}}\mathrm{\Delta }x}}\right)\left(F_{x_{i-\mathrm{1}/\mathrm{2},j,k}}^t-F_{x_{i+\mathrm{1}/\mathrm{2},j,k}}^t\right){\displaystyle }& {\displaystyle }+\left({\displaystyle \frac{\mathit{\omega }\mathrm{\Delta }t}{c_{n_{i,j,k}}\mathrm{\Delta }y}}\right)\left(F_{y_{i,j-\mathrm{1}/\mathrm{2},k}}^t-F_{y_{i,j+\mathrm{1}/\mathrm{2},k}}^t\right)+\left({\displaystyle \frac{\mathit{\omega }\mathrm{\Delta }t}{c_{n_{i,j,k}}\mathrm{\Delta }z}}\right)\\ {\displaystyle }& {\displaystyle }\left(F_{z_{i,j,k-\mathrm{1}/\mathrm{2}}}^t-F_{z_{i,j,k+\mathrm{1}/\mathrm{2}}}^t\right).\end{array}$$

The coefficients of Eq. (24) described in Eqs. (25)–(32) are for an internal grid cell with six neighbors. The coefficients for the top grid cells are modified for the presence of snow and/or standing water. A no-flux boundary condition was applied on the bottom grid cells; thus no geothermal flux was accounted for in this study. The coefficients for the grid cells on the lateral boundary are modified for a no-flux boundary condition. ELM 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 et al. (2013) for details about the computation of thermal properties and phase transition.

2.3.3 PETSc numerical solution

ELMv0, which considers flow only in the vertical direction, solves a tridiagonal and banded tridiagonal system of equations for water and energy transport, respectively. In ELM-3D, accounting for lateral flow in the subsurface results in a sparse linear system, Eqs. (10) and (24), where the sparsity 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. A description of the numerical tests that were conducted to ensure that the lateral coupling of hydrologic and thermal processes was correctly implemented is presented in the Supplement (Figs. S1 and S2)

2.4 Snow model and redistribution

The snow model in ELM-3D is the same as that in the default ELMv0 and CLM4.5 (Anderson, 1976; Dai and Zeng, 1997; Jordan, 1991), except for the inclusion of snow redistribution (SR). The snow model allows for a dynamic snow depth and up to five 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 et al., 2013). ELMv0 uses the methodology of Swenson and Lawrence (2012) to compute fractional snow cover area, which is appropriate for ESM-scale grid cells (∼ 100 km × 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 SR include topography and surface wind (Warscher et al., 2013); previous SR models include mechanistically based (Bartelt and Lehning, 2002; Liston and Elder, 2006) and empirically based (Frey and Holzmann, 2015; Helfricht et al., 2012) approaches. To mimic the effects of wind, we used a conceptual model to simulate 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.

2.5 System characterization

Hydrologic and thermal properties differ by depth and landscape type. We used the horizontal distribution of organic matter content from Wainwright et al. (2015) to infer soil hydrologic and thermal properties following the default representations in ELM. 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 for a given soil layer when its temperature falls below a threshold (T_threshold) of −2 ^∘C. This default value leads to overly large predicted latent and sensible heat fluxes during winter, compared to nearby eddy covariance measurements. We modified T_threshold to be 0 ^∘C in this study, resulting in improved predicted wintertime latent heat fluxes compared to the default version of the model (Fig. S3). 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.

2.6 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 (Fig. 1). The transect was 104 m long and 45 m deep and was discretized horizontally with a grid spacing of 0.25 m and an exponentially varying layer thickness in the vertical with 30 soil layers. The transect does not align with the sensor locations because our objective was not to validate the model for a few grid cells, but to focus on relative differences between predictions for rims and centers of a polygon field. No flow conditions for mass and energy were imposed on the east, west, and bottom boundaries of the domain. Temporal discretization of 30 min was used in the simulations. All simulations were performed in the satellite phenology (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 (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 (Fig. 1), and they were combined to analyze comparable landscape positions in the simulations (Fig. 2).

Figure 2 Simulated average winter snow surface elevation across the transect for three scenarios: (1) snow redistribution (SR) turned off and 1-D subsurface physics, (2) snow redistribution turned on and 1-D subsurface physics, and (3) snow redistribution turned on and 2-D subsurface physics. Surface elevation of the transect is shown by the solid black line. The dashed line indicates the boundary for comparison to observations in relatively lower (centers) and relatively higher (rims) topographical positions.

After testing, the model was used to investigate the effects of snow redistribution and 2-D subsurface hydrologic and thermal physics by analyzing three scenarios: (1) no snow redistribution and 1-D physics; (2) snow redistribution and 1-D physics; and (3) snow redistribution and 2-D 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.

3 Results and discussion

3.1 Snow depth

In the absence of SR, predicted snow depth exactly follows the topography. With SR, a much smaller dependence of winter average snow depth on topography is predicted (Fig. 2). Further, for the winter average, there are very small differences in snow depth between simulations with SR and 1-D or 2-D subsurface physics representations. Compared to observations, considering SR 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 (Fig. 3). There was no discernible difference in snow depth bias between the 1-D and 2-D physics (Table 1), although the predicted subsurface temperature fields were different, as shown below.

Table 1 Bias, root mean square error (RMSE), and correlation (R²) between modeled and observed snow depth at polygon center and rim, and the difference between center and rim for 2013 for three cases: snow redistribution (SR) off and 1-D physics, SR on and 1-D physics, and SR on and 2-D physics.

Figure 3 Monthly mean comparison of observation and simulated snow depth (a) in polygon rim and (b) in polygon center; (c) difference between polygon center and rim for 2013.

The temporal variation of the mean snow depth (Fig. 4a) and its spatial standard deviation (Fig. 4b) also differed based on whether SR was considered, but it was not affected by considering 2-D thermal or hydrologic physics. With SR, the snow depth coefficient of variation (Fig. 4c) was about 0.5 from December through the beginning of the snowmelt period, indicating relatively large spatial heterogeneity. Simulated snow depths for the three simulation scenarios are included in the Supplement (4).

Figure 4 (a) Mean, (b) standard deviation, and (c) coefficient of variation of simulated snow depth across the entire domain for 1-D and 2-D subsurface physics.

3.2 Soil temperature and active layer depth

Broadly, ELM-3D accurately predicted the polygon center soil temperature at depth intervals corresponding to the temperature probes (0–20, 20–50, 50–75, and 75–100 cm; Fig. 5a). Recall that the observed temperatures for the polygon center and rims were taken at single points in site A (Fig. 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 snowpack. 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 ELM-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 the Supplement (Figs. S5–S8).

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

Similarly, the soil temperatures were accurately predicted in the polygon rims (Fig. 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 snowpack acts as a heat insulator and keeps soil temperature higher in winter as compared to the polygon rims.

Three recent studies have used other mechanistic models to simulate soil temperature fields at this site and achieved comparably good comparisons with observations (Kumar et al., 2016 applied a 3-D version of PFLOTRAN; Atchley et al., 2015 and Harp et al., 2016 applied a 1-D version of the Arctic Terrestrial Simulator – ATS). However, those models used measured soil temperatures near the surface as the top boundary condition. In contrast, the top boundary condition in this work is the climate forcing (air temperature, wind, solar radiation, humidity, precipitation), and the ground heat flux is predicted based on ELM's vegetation and surface energy dynamics. We note that no parameter calibration was done in this work or that of Kumar et al. (2016), while the ATS parameterizations were calibrated to match the soil temperature profile.

Snow redistribution impacts spatial variability of soil temperature throughout the soil column. Absence of SR results in no significant spatial variability of soil temperature (Fig. 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 (Fig. 6b). With SR, energy dissipation in the lateral direction reduces the penetration depth of the soil temperature spatial variance (compare Fig. 6c and b).

Table 2 Bias, root mean square error (RMSE), and correlation (R²) between modeled and observed soil temperature at polygon center and rim at multiple soil depth for 2013 for three cases: snow redistribution (SR) off and 1-D physics, SR on and 1-D physics, and SR on and 2-D physics.

Figure 6 Simulated daily spatial standard deviation for each soil layer averaged across 10 years of near-surface soil temperature for the simulation performed with snow redistribution turned off and 1-D subsurface physics (a); snow redistribution turned on and 1-D subsurface physics (b); and snow redistribution turned on and 2-D subsurface physics (c).

With 1-D physics, the average spatial and temporal difference of the active layer depth (ALD) between simulations with and without SR was 1.7 cm (Fig. 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 1-D version of ATS that SR resulted in deeper thaw depths in both polygon centers (by ∼ 3 cm) and rims (∼ 0.3 cm). Thus, their results for polygon centers are consistent in sign but lower in magnitude than ours, but opposite in sign for the rims.

Figure 7 Temporal mean of the bottom of the active layer (a) and standard deviation of the active layer depth (b) over the 10-year period across the modeling domain.

Across 10 years of simulation, the interannual variability (IAV) in ALD varied substantially between the three scenarios (Fig. 7b). As expected, for the 1-D physics without SR scenario (green line), the IAV in ALD was determined by a landscape position because of differences in soil and vegetation parameters. With SR and 1-D 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 1-D versus 2-D physics on the ALD across the transect was modest (mean absolute difference ∼ 3 cm). Generally, because 2-D 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; Fig. 7a) than with 1-D physics. This difference was also reflected in the thaw depth IAV across the transect, where 2-D physics led to a smoother lateral profile of interannual variability than with 1-D physics.

The impact of physics formulation (i.e., 1-D or 2-D) alone was investigated by analyzing differences between soil temperature profiles over time for polygon rims and centers in simulations with snow redistribution. Inclusion of 2-D subsurface physics resulted in soil temperatures with depth and time that were lower in the polygon rims (Fig. 8a) and higher in polygon centers (Fig. 8b). Using the simulations from the scenario with SR and 2-D physics, we evaluated the extent to which 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 (Fig. 9). During winter, both polygon rims and centers were predicted to have soil temperature spatial variability > 1 ^∘C up to a depth of ∼ 2 m. The soil temperature spatial variability in winter due to snow redistribution was dissipated over the summer. During the summer, polygon centers were relatively more homogeneous vertically compared to polygon rims.

Figure 8 Time series of spatial mean soil temperature differences between “SR = on + physics = 1-D” and “SR = on + physics = 2-D” at polygon rim (a) and polygon center (b).

3.3 Surface energy budget

Predicted monthly mean and spatial mean (μ) surface latent heat fluxes across the transect were very similar between the three scenarios (Fig. 10a), with a growing seasonal mean difference of < 1.0 W m⁻². However, the spatial variability (SV =σ; Fig. 10b) and coefficient of variation (CV $=\mathit{\sigma }/\mathit{\mu }$; Fig. 10c) of latent heat fluxes were different between the scenarios with SR (1-D and 2-D 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.

Figure 9 Time series of soil temperature spatial standard deviation for “SR = on + physics = 2-D” at polygon rim (a) and polygon center (b).

Figure 10 Latent heat flux interannual (a) mean, (b) standard deviation, and (c) coefficient of variation across the site A transect.

The predicted temporal monthly mean and spatial mean surface sensible heat fluxes across the transect were also similar between the three scenarios (Fig. 11a), with a growing season mean absolute difference of < 3.5 W m⁻². 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⁻² (Fig. 11b) and CVs of ∼ 1.5 (Fig. 11c). As for the latent heat fluxes, the differences in standard deviation and CV of sensible heat fluxes were small between the 1-D and 2-D scenarios with SR, arguing that the subsurface lateral energy exchanges associated with the 2-D 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⁻² in May; Fig. 10b).

Figure 11 Same as Fig. 10 except for sensible heat flux.

3.4 Soil moisture

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

3.5 Caveats and future work

The good agreement between ELM-3D predictions and soil temperature observations demonstrates 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.

Figure 12 Same as Fig. 6 except for liquid saturation.

ELMv0, a one-dimensional model, is embarrassingly parallel with no cross-processor communication. The current implementation of the three-dimensional solver in ELM-3D only supports serial computing. Support of parallel computing will be included in a future version of the model. Because of computational constraints, we applied a 2-D transect domain to the site, instead of a full 3-D domain. We are working to improve the computational efficiency of the model, which will facilitate a thorough analysis of the effects of 3-D 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 2-D surface thermal flow model in a 3-D domain, which is another goal for our future work. However, we note that the good agreement using the 2-D 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 ELM, although incorporation of these processes is a long-term development goal.

The current representation of vegetation in ELM-3D for these polygonal tundra systems is oversimplified. 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 on each plant functional type in the domain, ignores the likely influence of nutrient constraints (Zhu et al., 2016) 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 et al., 2015).

Development of subgrid parameterizations to parsimoniously capture fine-scale processes will be pursued in the future. For example, a two-tile approach to represent hydrologic and thermal processes in coupled polygon rims and centers with snow redistribution should be evaluated. Inclusion of lateral subsurface processes has a greater impact on predicted subgrid variability than on spatially averaged states. Thus, one possible extension of the current model would be to explicitly include an equation for the temporal evolution of subgrid variability using the approach of Montaldo and Albertson (2003). The use of reduced-order models (e.g., Pau et al., 2014; Liu et al., 2017) is an alternate approach to estimate fine-scale hydrologic and thermal states from a coarse resolution representation. Additionally, lateral subsurface processes can be included in the land surface model via a range of numerical discretization approaches of varying complexity, e.g., adding lateral water and energy fluxes as source–sink terms in the existing 1-D model, implementing an operator split approach to solve vertical and lateral processes in a non-iterative approach, or solving a fully coupled 3-D model. Trade-offs between approaches that represent lateral processes and computational costs need to be carefully studied before developing quasi or fully three-dimensional land surface models. While the present study focused on application and validation of ELM-3D at fine-scale, future work will focus on regional-scale applications using comprehensive datasets and the Distributed Model Intercomparison Project Phase 2 modeling protocol (Smith et al., 2012). Although we found no significant effect of topography and SR on the 100 m × 100 m grid-averaged exchanges with the atmosphere, future work needs to analyze intermediate-scale (e.g., 100 m–10 km) topographical variation and the potential effects on biogeochemical and plant processes and surface exchanges.

4 Summary and Conclusions

In a polygonal tundra landscape, we analyzed effects of microtopographical surface heterogeneity and lateral subsurface transport on soil temperature, soil moisture, and surface energy exchanges. Starting from the climate-scale land model ELMv0, we incorporated in ELM-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.

Over the observational record, ELM-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² of 0.59 ^∘C, 1.82 ^∘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 maximum thaw depth under the polygon centers. The effect of lateral energy fluxes on active layer depths was ∼ 3 cm. Compared to 1-D physics, the 2-D 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 1-D columns. During the summer, polygon center temperatures were relatively more vertically homogeneous than temperatures in the rims.

The monthly mean 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 overprediction 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.

Overall, our analysis demonstrates the potential and value of explicitly representing snow redistribution and lateral subsurface hydrologic and thermal dynamics in polygonal ground systems and quantifies the effects of these processes on the resulting system states and surface energy exchanges with the atmosphere. The integration of a 3-D subsurface model in the E3SM Land Model also allows for a wide range of analyses heretofore impossible in an Earth System Model context.