the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
New insights into the drainage of inundated icewedge polygons using fundamental hydrologic principles
Vitaly Zlotnik
Charles J. Abolt
Bob Busey
Sofia T. Avendaño
Brent D. Newman
Adam L. Atchley
Elchin Jafarov
Cathy J. Wilson
Katrina E. Bennett
The pathways and timing of drainage from the inundated centers of icewedge polygons in a warming climate have important implications for carbon flushing, advective heat transport, and transitions from methane to carbon dioxide dominated emissions. Here, we expand on previous research using a recently developed analytical model of drainage from a lowcentered polygon. Specifically, we perform (1) a calibration to field data identifying necessary model refinements and (2) a rigorous model sensitivity analysis that expands on previously published indications of polygon drainage characteristics. This research provides intuition on inundated polygon drainage by presenting the first indepth analysis of drainage within a polygon based on hydrogeological first principles. We verify a recently developed analytical solution of polygon drainage through a calibration to a season of field measurements. Due to the parsimony of the model, providing the potential that it could fail, we identify the minimum necessary refinements that allow the model to match water levels measured in a lowcentered polygon. We find that (1) the measured precipitation must be increased by a factor of around 2.2, and (2) the vertical soil hydraulic conductivity must decrease with increasing thaw depth. Model refinement (1) accounts for runoff from rims into the icewedge polygon pond during precipitation events and possible rain gauge undercatch, while refinement (2) accounts for the decreasing permeability of deeper soil layers. The calibration to field measurements supports the validity of the model, indicating that it is able to represent icewedge polygon drainage dynamics. We then use the analytical solution in nondimensional form to provide a baseline for the effects of polygon aspect ratios (radius to thaw depth) and coefficient of hydraulic conductivity anisotropy (horizontal to vertical hydraulic conductivity) on drainage pathways and temporal depletion of ponded water from inundated icewedge polygon centers. By varying the polygon aspect ratio, we evaluate the relative effect of polygon size (width), interannual increases in activelayer thickness, and seasonal increases in thaw depth on drainage. The results of our sensitivity analysis rigorously confirm a previous analysis indicating that most drainage through the active layer occurs along an annular region of the polygon center near the rims. This has important implications for transport of nutrients (such as dissolved organic carbon) and advection of heat towards icewedge tops. We also provide a comprehensive investigation of the effect of polygon aspect ratio and anisotropy on drainage timing and patterns, expanding on previously published research. Our results indicate that polygons with large aspect ratios and high anisotropy will have the most distributed drainage, while polygons with large aspect ratios and low anisotropy will have their drainage most focused near their periphery and will drain most slowly. Polygons with small aspect ratios and high anisotropy will drain most quickly. These results, based on parametric investigation of idealized scenarios, provide a baseline for further research considering the geometric and hydraulic complexities of icewedge polygons.
Polygonal tundra occurs in continuouspermafrost landscapes lacking exposed bedrock or active sedimentation (Brown et al., 1997; Mackay, 2000). Estimates of its spatial extent vary from around 250 000 km^{2} (Minke et al., 2007), or approximately the size of England, to 2.6 million km^{2}, or approximately 30 % of the Arctic land surface (Mackay, 1972). This terrain is rich in organic carbon (Tarnocai et al., 2009; Hugelius et al., 2014), which is poised for release to the atmosphere as carbon dioxide or methane under a warming climate and may cause a significant feedback to climate change (Schuur et al., 2008). Depressions in the microtopography of polygonal tundra collect water from precipitation and snowmelt events, resulting in a slow release from the landscape through the thawed subsurface to surface water drainage networks (Helbig et al., 2013). This process has important implications for the leaching of dissolved organic carbon from activelayer soils and for advective heat transfer through the subsurface. The role that polygonal tundra regions play in global terrestrial biogeochemical feedbacks as the Arctic warms necessitates a greater understanding of their hydrologicflow regimes.
Polygonal tundra forms in cold environments by the cyclic process of vertical cracking of frozen ground due to thermal contraction, water infiltration into these cracks, freezing of this infiltrated water, and subsequent recracking. Over many cycles, this process leads to the growth of subsurface ice wedges connected in polygonal patterns known as icewedge polygons (Liljedahl et al., 2016). Icewedge polygons vary in size, from several meters to a few tens of meters in diameter, and develop over time frames of hundreds to thousands of years (Leffingwell, 1915; Lachenbruch, 1962; Mackay, 2000; Abolt and Young, 2020). Thermal expansion of the ice and soil during summer months often results in warping of the soil strata, forming parallel rims on both sides of the ice wedge and a trough directly above the ice wedge (refer to Figs. 1 and 2). Icewedge polygons with wellformed rims and troughs with a distinct, central topographic depression are referred to as lowcentered polygons (as opposed to highcentered or transitional polygons). Ponding occurs frequently in the centers because the rims can trap water during snowmelt or precipitation events. This ponded water subsequently drains through the subsurface to the troughs (Helbig et al., 2013; Koch et al., 2018; Wales et al., 2020).
Recent observational studies have shaped our current conceptualization of lowcentered polygonal tundra hydrology (Boike et al., 2008; Helbig et al., 2013; Koch et al., 2014; Liljedahl et al., 2016; Koch, 2016; Koch et al., 2018), indicating a system of inundated polygonal centers surrounded by elevated rims that prevent surface water from flowing overland into surrounding trough ponds. This land formation results in standing water over the polygon centers, reduces immediate landscape runoff from precipitation events, and increases evaporation (Liljedahl et al., 2016). Additionally, the elevated ponded water in polygon centers produces hydrological gradients that result in sustained outward flow through the subsurface under the rims, a process that has been observed at several field sites (Helbig et al., 2013; Liljedahl and Wilson, 2016; Koch, 2016; Wales et al., 2020). This lateral flow controls landscape redistribution of water during the summer months (Helbig et al., 2013), governs ponded water budgets (Koch et al., 2014; Koch, 2016), and makes up a notable portion of regional river discharge (King et al., 2020). The manner and timing with which polygonal tundra landscapes transition from inundated to drained conditions have important implications for (1) transitions from atmospheric emissions of methane to carbon dioxide (Conway and Steele, 1989; Moore and Dalva, 1997; Zona et al., 2011; Zhu et al., 2013; O'Shea et al., 2014; Throckmorton et al., 2015; Wainwright et al., 2015; Lara et al., 2015; Vaughn et al., 2016), (2) dissolved organic carbon emissions to surface waters (Zona et al., 2011; Abnizova et al., 2012; Laurion and Mladenov, 2013; Larouche et al., 2015; Plaza et al., 2019), (3) biological succession (Billings and Peterson, 1980; Jorgenson et al., 2015; Wolter et al., 2016), and (4) ground surface deformation (Mackay, 1990, 2000; Raynolds et al., 2014; Nitzbon et al., 2019).
The microtopographic features of polygonal tundra result in pronounced finescale spatial gradients in thermal and hydrologic conditions (Boike et al., 2008; Zona et al., 2011; Helbig et al., 2013), leading to sharply contrasting biogeochemical conditions (Newman et al., 2015; Norby et al., 2019). Standing water in polygonal centers shows distinctly different chemistry compared to surface water in troughs and ponds (Zona et al., 2011; Koch et al., 2014, 2018), which has been attributed to shallow subsurface flow from centers to troughs (Koch, 2016). For example, Koch et al. (2014) observed high nutrient concentrations in troughs adjacent to centers with low concentrations, indicating drainage from the latter to the former. Given that the icewedge polygon can be considered to be the fundamental hydrologic landscape unit that initiates landscapescale surface or subsurface flow and discharge, understanding their characteristic drainage pathways and residence times provides a basis to quantify the transport of terrestrially sourced dissolved organic carbon to surface waters. This is important for carbon emissions because dissolved organic carbon in surface waters can be mineralized and released to the atmosphere (Raymond et al., 2013).
Although there have been numerous field observations of inundated lowcentered polygon drainage to troughs (Helbig et al., 2013; Koch et al., 2014; Wales et al., 2020), the research here is, to our knowledge, the first indepth investigation into the relative pathways and timing associated with inundated icewedge polygon drainage based on hydrologic first principles. Aside from the field tracer experiments of Wales et al. (2020), previous research has primarily focused on 1D vertical hydrothermal effects within icewedge polygons (Atchley et al., 2015; Harp et al., 2016; Atchley et al., 2016) or polygonal tundra surface hydrology (Boike et al., 2008; Jan et al., 2018a, b). Other researchers have investigated the hydrology of multiple polygons without investigating the drainage pathways within individual icewedge polygons (CrestoAleina et al., 2013; Nitzbon et al., 2019, 2020). While basic hydrology dictates that icewedge polygon geometry and heterogeneity will explicitly govern subsurface drainage pathways and duration, to date, a rigorous hydrological investigation of this process has not been presented.
The geometric shape of lowcentered polygons along with soil hydraulic properties (for example, hydraulic conductivity) affects the distribution of hydraulic heads that control the pathways and timing of inundated icewedge polygon drainage. In this paper, we use a recently developed model (Zlotnik et al., 2020) based on hydrogeological first principles (Harr, 1962; Cedergren, 1968; Freeze and Cherry, 1979; Bear, 1979) to understand and gain intuition into the effects of geometry and hydraulic properties on inundated icewedge polygon drainage.
This model idealizes the thawed subsurface of an inundated lowcentered polygon as a thin cylinder overlain with an initial height of ponded water that drains through the cylindrical thawed soil layer to the surrounding trough (outer vertical boundary; an annular ring defined by a line from (R,0) to (R,L) in Fig. 1). The depth (L) and radius (R) of the thawed soil layer of the polygon center and horizontal (radial) and vertical hydraulic conductivity are parameters of the solution. It is important to note that we are focusing on the drainage of inundated lowcentered polygons here and that biogeochemical investigations suggest that the drainage of noninundated lowcentered, transitional, and highcentered polygons may have different characteristics (Heikoop et al., 2015; Newman et al., 2015; Wainwright et al., 2015; Wales et al., 2020).
We demonstrate that the model is able to accurately simulate icewedge polygon drainage by performing a calibration to field measurements over an entire thaw season. We calibrate the model parameters to fit water levels measured in the center of a lowcentered polygon within the Barrow Environmental Observatory near Utqiaġvik, Alaska. We force the model with measured trough water levels, precipitation, thaw depth, and evapotranspiration. The goodness of fit of initial calibration efforts indicated that the model was unable to capture the polygoncenter water level trends (i.e., the model was falsified). However, by considering (1) increased precipitation and (2) a thawdepthdependent (and therefore timevariable) vertical hydraulic conductivity, the model produces a good fit to the water levels. Model refinement (1) accounts for pond accumulation due to runoff from surrounding microtopographic highs (e.g., rims), while model refinement (2) accounts for the decreasing hydraulic conductivity of deeper tundra soil layers. The calibration indicates the minimum refinements necessary to a hydrology model to capture icewedge polygon drainage dynamics and indicates the important factors in this process. This calibration, based on transient boundary conditions, provides confidence in the physical meaningfulness of the steadystate snapshots presented in the sensitivity analysis.
We then use the model to investigate the pathways and timing of inundated polygon drainage through a sensitivity analysis of polygon geometry and anisotropy. The analysis of the model enhances our intuition into the pathways and timing of inundated icewedge polygon drainage for polygons with different geometries and degrees of anisotropy. Varying the geometry of the cylinder allows us to capture not only relative differences in drainage between differentsized polygons (in other words, polygons with different radii) but also seasonal and interannual variations in polygon thawlayer depths within a single polygon (in other words, as the thawed soil layer increases during a single thaw season or as the active layer increases from year to year). By allowing for different effective hydraulic conductivities between the horizontal and vertical directions, we can evaluate the effects of preferential flow directions on drainage. In addition to geologic layering, processes such as frost heave and horizontal ice lens thaw can result in preferential horizontal flow in cold climates (Mackay, 1981; Matsuoka and Moriwaki, 1992; Wales et al., 2020). It has also been observed that highhydraulicconductivity peat layers overlying lowerhydraulicconductivity mineral soil result in horizontal watershed drainage (McDonnell et al., 1991; Brown et al., 1999; Quinton and Marsh, 1999). Helbig et al. (2013) found that in polygonal tundra, the peat layer quickly redistributed precipitation events across the subsurface topography from the rims to the centers and troughs. Vertical cracks, voids left from decayed roots, and animal burrows can result in preferential vertical flow. The cylindrical geometry and anisotropy are not intended to cover all potential variations in icewedge polygons but rather provide base cases where those variations will cause deviations around the idealizations considered here.
The calibration to field data indicates that the model is able to capture icewedge polygon drainage dynamics. Our sensitivity analysis provides a new perspective on inundated polygon hydrogeology, indicating the effects of polygon geometry and hydraulic conductivity anisotropy on drainage pathways and timing. Although the simplifications of the model may limit its applicability to some scenarios, they allow general intuitive insights to be drawn which would be obfuscated without them. The findings here provide a basis to quantify and understand deviations from our idealized scenarios.
2.1 Model overview
We use recently developed 3D axisymmetric analytical solutions, derived and validated in Zlotnik et al. (2020) and described in Appendix A, of hydraulic heads and the stream function in the thawed soil layer below a polygon center (dashed rectangle in Fig. 1) to investigate inundated lowcentered polygon drainage pathways. While we perform the calibration to field data using the dimensional form of the solution for hydraulic heads, we perform the sensitivity analyses using the nondimensional form of the solution for hydraulic heads and stream function (refer to Eqs. A7 and A8). The nondimensional hydraulic heads and stream function are independent of time, representing the relative pattern of hydraulic heads and stream function throughout the drainage process. The dimensional values at different times would indicate the change in the absolute magnitude of these patterns as the drainage process proceeds. To evaluate drainage pathways, we plot flow nets, as illustrated in Fig. 1, composed of lines of equal nondimensional hydraulic heads (${h}^{*}({r}^{*},{z}^{*})$), referred to as equipotentials, and contours of the stream function (${\mathrm{\Psi}}^{*}({r}^{*},{z}^{*})$) as a function of radius and depth. The stream function contours are referred to as streamlines and are used here to identify steady drainage pathways.
Nondimensional hydraulic head lines are drawn from 0.05 to 0.95 with increments of 0.1. Dimensional heads ($h(r,z,t)$) can be obtained from the nondimensional heads (${h}^{*}({r}^{*},{z}^{*})$) using the ponded height of water in the polygon center H_{c}(t) and the height of water in the trough H_{t} as
where r^{*} and z^{*} are nondimensional radius and depth, respectively; L is the polygon thaw depth; and K_{r} and K_{z} are the radial (horizontal) and vertical hydraulic conductivity, respectively (refer to Fig. 1). The aspect ratio is defined as $R/L$, and coefficient of anisotropy is ${K}_{\mathrm{r}}/{K}_{\mathrm{z}}$.
We use a Robin boundary condition (a thirdtype boundary condition allowing both head and flux to be specified) for the outer vertical boundary. This allows the model to represent the resistance of drainage under the rim due to elevated frozen ground and due to the accumulation of fine soil particles at the soil–water interface of the trough (Koch et al., 2018). The resistance to drainage below the polygon rim and across the soil–water interface into the trough is represented by a hydraulic conductance $\mathit{\kappa}=k/l$, where k is the hydraulic conductivity, and l is the thickness of the drainage interface (refer to Eq. A3).
We normalize the stream function $\mathrm{\Psi}({r}^{*},{z}^{*})$ from 0 to 1 and, similar to heads, plot stream function contours (streamlines) from 0.05 to 0.95 with increments of 0.1. The region between any two adjacent streamlines conveys 10 % of the drainage, while the two remaining regions (less than 0.05 and greater than 0.95) each convey 5 % of the drainage. To make quantitative comparisons of the relative spread of drainage between polygons, we calculate the percent of each polygon volume accessed by 95 % of the drainage (polygon volume with nondimensional stream function >0.05).
The change in nondimensional ponded water height in the polygon center due solely to drainage (the depletion curve) over time can be expressed by a simple exponential decay function as
where t_{L} is the characteristic time of drainage (refer to Eq. A11), defined as the time when the ponded height is $\mathrm{1}/e$ times, or ∼ 37 % of, its original height H_{c,0}. The dimensional ponded height H_{c}(t) can be obtained as
As described in Zlotnik et al. (2020), for the calibration, temporally changing trough water levels, precipitation, evapotranspiration, and thaw depth are accounted for by running the model in a piecewise fashion.
2.2 Calibration approach
We performed nonlinear leastsquares calibrations using a Levenberg–Marquardt approach (Transtrum et al., 2011). The objective function is the sumofsquared errors (SSEs) between the measured and simulated polygoncenter water levels expressed as
where θ is a vector of model parameters, and ${H}_{\mathrm{c}}^{m}\left({t}_{i}\right)$ is the measured polygoncenter water level at the ith time. We present three calibration cases, including (1) a base case (model setup similar to the validation in Zlotnik et al., 2020), (2) precipitation multiplied by a calibrated factor, and (3) a depthdependent vertical hydraulic conductivity (K_{z}=f(D)). Depending on the calibration case, θ contains different parameters.
In calibration case 1, θ includes parameters K_{r}, K_{z}, κ, and H_{c,0}. Based on field observations and the validation in Zlotnik et al. (2020), we constrain K_{z} to be less than K_{r} using a metaparameter ${\mathcal{F}}_{{K}_{\mathrm{z}}}$ as ${K}_{\mathrm{z}}={K}_{\mathrm{r}}\mathit{\sigma}\left({\mathcal{F}}_{{K}_{\mathrm{z}}}\right)$, where σ is the sigmoid function defined as $\mathit{\sigma}\left(x\right)=\mathrm{1}/(\mathrm{1}+\mathrm{exp}(x\left)\right)$.
In calibration case 2, we add a precipitation multiplier M_{P} to case 1, where $\widehat{P}={M}_{\mathrm{P}}P$, and $\widehat{P}$ is the augmented precipitation accounting for runoff from microtopographic highs and potential rain gauge undercatch.
In calibration case 3, we add vertical hydraulic conductivity as a superelliptical function of thaw depth to calibration case 2 as ${K}_{\mathrm{z}}\left(D\right)={K}_{\mathrm{z}}^{\mathrm{min}}+({K}_{\mathrm{z}}^{\mathrm{max}}{K}_{\mathrm{z}}^{\mathrm{min}}){\left(\mathrm{1}{D}_{\mathrm{norm}}^{a}\right)}^{\mathrm{1}/a},\phantom{\rule{0.33em}{0ex}}{D}_{\mathrm{norm}}=\frac{D{D}^{\mathrm{min}}}{{D}^{\mathrm{max}}{D}^{\mathrm{min}}}$, where ${K}_{\mathrm{z}}^{\mathrm{max}}$ and ${K}_{\mathrm{z}}^{\mathrm{min}}$ are the maximum and minimum values of K_{z}, respectively; D^{max} and D^{min} are the maximum and minimum values of D during the thaw season, respectively; and a is the superelliptical shape parameter that controls curvature. ${K}_{\mathrm{z}}^{\mathrm{max}}$ is constrained to be positive using a metaparameter ${\mathcal{F}}_{{K}_{\mathrm{z}}^{\mathrm{max}}}$ and the softplus function as ${K}_{\mathrm{z}}^{\mathrm{max}}=\mathrm{log}(\mathrm{1}+\mathrm{exp}({\mathcal{F}}_{{K}_{\mathrm{z}}^{\mathrm{max}}}\left)\right)$. We ensure that ${K}_{\mathrm{z}}^{\mathrm{min}}\le {K}_{\mathrm{z}}^{\mathrm{max}}$ using a metaparameter ${\mathcal{F}}_{{K}_{\mathrm{z}}^{\mathrm{min}}}$ as ${K}_{\mathrm{z}}^{\mathrm{min}}={K}_{\mathrm{z}}^{\mathrm{max}}\mathit{\sigma}\left({\mathcal{F}}_{{K}_{\mathrm{z}}^{\mathrm{min}}}\right)$. We ensure that the shape parameter a is restricted between 0.5 and 1.5 using a metaparameter ℱ_{a} as $a=\mathrm{0.5}+\mathrm{1.5}\mathit{\sigma}\left({\mathcal{F}}_{a}\right)$. The constrained superelliptical function allows the calibration to identify a general nonlinear decrease in vertical hydraulic conductivity with increasing thaw depth.
2.3 Acquisition of field data used in calibration
The measured data used in the calibration were collected during the thaw season of 2013 from 16 June to 18 September from a lowcentered polygon within the Barrow Environmental Observatory near Utqiaġvik, Alaska (Fig. 2). We used water levels collected by Liljedahl and Wilson (2016) with 0.2 cm resolution and precipitation collected by Hinzman et al. (2014) with 0.1 mm resolution. We determined thaw depths from thermal transects collected by Romanovsky et al. (2017) with an accuracy of 0.1 ^{∘}C. Due to a lack of continuous local evapotranspiration measurements, we obtained evapotranspiration data from NASA’s Global Land Data Assimilation System (GLDAS) (Rodell et al., 2004). Given the continuous (albeit diurnally fluctuating) lowmagnitude evapotranspiration signal, its effect on our calibration is relatively insignificant compared to the sporadic precipitation events that drive largescale fluctuations in water levels. Therefore, in lieu of local evapotranspiration measurements, the GLDAS evapotranspiration is deemed sufficient for our purposes here.
2.4 Parameter ranges for sensitivity analyses
We selected aspect ratio and anisotropy scenarios based on existing literature and observations. While a panArctic survey of icewedge polygon diameters does not to our knowledge currently exist, researchers have provided general characteristics based on extensive observations. Leffingwell (1915) states that icewedge polygons have an estimated average diameter of around 15 m. Liljedahl et al. (2016) indicate that lowcentered polygons (including rims) can have diameters from 5–30 m. Abolt et al. (2018) state that polygonal formations in icewedge tundra are 10–30 m in diameter. Based on digital elevation models, Abolt et al. (2019) calculate that the mean polygon diameter near Utqiaġvik, Alaska, is around 15 m, with 5th and 95th percentiles of approximately 8 and 33 m, respectively. Using a similar analysis on data from Prudhoe Bay (Abolt and Young, 2020), the mean polygon diameter is around 13 m, with 5th and 95th percentiles of approximately 7 and 30 m, respectively.
Maximum thaw depths are increasing throughout the Arctic (Jorgenson et al., 2006). Deeper thaw depths may ultimately result in lowcentered polygon transition to highcentered polygons, at which time polygoncenter inundation will cease to occur due to rim collapse. Therefore, there is a practical limit to the thaw depth of interest in our analysis. Lewkowicz (1994) reported that activelayer thickness varied from 40 to 90 cm within the polygonal tundra of Fosheim Peninsula, Ellesmere Island, Canada, in 1994. Shiklomanov et al. (2010) reported that polygonal tundra at several study sites near Utqiaġvik, Alaska, monitored from 1995–2009 had activelayer thicknesses from 17–47 cm. Based on extensive groundpenetrating radar surveys near Utqiaġvik, Alaska, in 2013, Jafarov et al. (2017) document that average activelayer thickness was 41 cm with a standard deviation of 9 cm. Even for polygons with substantial activelayer thickness, it is important to understand the change in drainage throughout the season by evaluating earlyseason thaw depths, which are captured in our analysis by larger aspect (i.e., radius to thaw depth) ratios.
Based on these considerations, we evaluate aspect ratios from 2.5 to 20. For example, given a thaw depth of 1 m, an aspect ratio of 2.5 would represent a 2.5 m radius polygon center, while an aspect ratio of 20 would represent a 20 m radius polygon center, more than covering the observed range. To consider cases with thinner thawed soil layers, larger aspect ratios can be used. For example, given a thaw depth of 0.5 m, an aspect ratio of 20 would represent a polygon with a 10 m radius center.
Comprehensive measurements of hydraulic conductivity anisotropy are lacking from icewedge polygons. Based on tracer arrival times, Wales et al. (2020) estimated horizontal conductivities of approximately 0.7–84 m d^{−1} for a lowcentered polygon and 0.01–0.3 m d^{−1} for a highcentered polygon. Wales et al. (2020) stated that these are considered lowerbound estimates of horizontal conductivity since, based solely on breakthrough times, the approach is unable to isolate horizontal and verticalflow effects on arrival times. Beckwith et al. (2003) performed laboratory measurements of anisotropy on small samples from peat soils where horizontal conductivities were around 6–7 m d^{−1}, and vertical conductivities were around 0.2–0.4 m d^{−1}, indicating preferential horizontal flow. The soil samples used in Beckwith et al. (2003) were from peat bogs in England and therefore would not have been subjected to similar freeze–thaw dynamics as contemporary polygonal tundra soils, which would presumably lead to even greater preferential horizontal flow. Therefore, these estimates may provide a lowerbound estimate for Arctic peat horizontal hydraulic conductivity. Based on the existing literature and above considerations, we considered hydraulic conductivities from 0.005–5 m d^{−1}, primarily focusing on anisotropic values from 0.1–100. These values allow us to investigate the relative effects of anisotropy on the drainage pathways and timing using hydraulic conductivities consistent with currently available measurements.
A physical interpretation of our selected values of anisotropy coefficient can be obtained by considering that icewedge polygon soils are typically layered and that the horizontal and vertical hydraulic conductivities can therefore be considered to be effective properties. As such, the effective horizontal hydraulic conductivity captures parallel flow dominated by the higherconductivity layers, and the effective vertical hydraulic conductivity captures flow in series across multiple layers and is dominated by the lowerconductivity layers. Therefore, an anisotropy coefficient of 100 would effectively capture layers with 2 orders of magnitude difference in hydraulic conductivity, while an anisotropy coefficient of 0.1 would capture the hypothetical scenario where vertical cracks or burrows result in preferential vertical flow. Given the current lack of direct measurements of icewedge polygon anisotropy coefficient, we cover a broad range of possible scenarios.
The hydraulic conductivity of the interface between soil and open water can be half that of the rest of the soil due to the accumulation of fines (Koch et al., 2018). The hydraulic conductance under the rim may also impede drainage due to a raised permafrost table following the surface topography. Note that while the raised permafrost table under the rim will constrict flow and alter drainage pathways, hydrologic first principles indicate that this effect will be restricted to the region of the model near the rim. This is analogous to the effect of a partially penetrating well on flow in an aquifer, which dissipates quickly and is nonexistent by a lateral distance of 1.5 to 2 times the aquifer thickness for isotropic aquifers (Bear, 1979). This effect will diminish with increasing aspect ratio and decreasing anisotropy coefficient and will not alter the qualitative insights drawn from the relative comparison of drainage pathways in our analysis. It is conceivable that early in the thaw season the permafrost table under the rim could extend above the center ground surface within the vertical interval of the center pond. In this case, there will likely be very little drainage occurring associated with a very small discharge conductance. Although the conductance of the outer interface (defined as κ above) can therefore influence drainage, a preliminary analysis indicated that in most cases, its effect on drainage was significantly less than that of aspect ratio and anisotropy coefficient (Zlotnik et al., 2020). Therefore, we chose to use a default value for κ in all cases (unless otherwise specified) of 5 d^{−1}.
In practice, as in our calibration, the water level in troughs (H_{t} in Eqs. 1 and 3) will vary throughout the thaw season, affecting the magnitude of heads in the soil of the polygon center and drainage times. As the nondimensional heads are relative to H_{t} (refer to Eq. A7), its value does not affect our relative comparisons of drainage patterns (which are based on nondimensional heads that are normalized from 0 to 1). The value of H_{t} will affect our comparisons of drainage time but in a systematic, interpretable manner. For example, a higher H_{t} will compress the exponential curve defined by Eq. (3) upwards, while a lower H_{t} will stretch the exponential curve downwards. In cases where H_{t} is below the polygoncenter ground surface, the solution is valid until H_{c} reaches the ground surface, at which time the ponded center has completely drained. Therefore, to isolate our analysis to aspect ratio and anisotropy, we have set H_{t} in all cases equal to the polygoncenter ground surface.
Although included in the calibration, we have neglected the effects of evaporation and precipitation in the sensitivity analysis as they will not affect the drainage patterns we present (based on nondimensional heads), and their effect on drainage timing (based on nondimensional depletion curves) is straightforward, shifting the nondimensional exponential drainage curve upwards or downwards. In other words, using nondimensional variables is a powerful approach to gain intuition into the fundamentals of inundated icewedge polygon drainage irrespective of variable magnitude.
We verify the model through calibration to water level measurements, identifying refinements necessary for hydrologic models to match field observations of polygon drainage. We present drainage flow nets for various aspect ratios and anisotropy coefficient values, polygoncenter ponded height depletion curves, and maps of the percent of the thawed soil accessed by 95 % of the drainage flow and depletion characteristic times as a function of aspect ratio and anisotropy coefficient. We compare the effect of aspect ratio on drainage pathways in an isotropic and a highly anisotropic (${K}_{\mathrm{r}}/{K}_{\mathrm{z}}=\mathrm{100}$) polygon. The drainage pathways are also evaluated for various anisotropies holding the aspect ratio constant. We present a global perspective of the combined effects of aspect ratio and anisotropy coefficient on the focusing and spreading of the drainage flow by mapping the percent of the polygon thawed soil volume that is accessed by 95 % of the drainage flow. We evaluate the effect of aspect ratio and anisotropy coefficient on drainage time by plotting the polygonponded water height depletion curves. Similar to accessed volume, we provide a global perspective on the combined effects of aspect ratio and anisotropy coefficient on drainage time by mapping the depletion characteristic time. We illustrate the counteracting effects of aspect ratio and anisotropy coefficient on drainage pathways by showing two, although geometrically and hydrologically dissimilar, mathematically equivalent solutions of drainage.
3.1 Calibration to field measurements
We present the bestfit polygoncenter water levels for the calibration cases along with the measured values in Fig. 3. The two dominant drivers of the polygoncenter water levels, the measured trough water level and precipitation, are shown for reference. Calibration case 1 (green line) is unable to match the high and low points of the water level measurements and instead follows a more medial path over time (RMSE = 1.42 cm, R^{2}=0.49; refer to Table 1). This indicates that the measured precipitation alone is not able to account for the increased water levels after rain events. We also consider the calibration to be failed (the model to be falsified; Popper, 2005) because the calibrated value of the discharge conductance (κ=835.3 d^{−1}; refer to Table 1) is extremely large and nonphysical.
Considering that precipitation will run off from rims and collect in the polygoncenter pond and that rain gauges may have undercatch issues (e.g., Pollock et al., 2018), we performed calibration case 2, also shown in Fig. 3 (purple line), which includes a precipitation multiplier. Calibration case 2 significantly improves the fit (RMSE = 0.89 cm, R^{2}=0.80) but still has significant mismatch as the limitations of the model require the calibration to compromise between matching water levels at early and late times. As in calibration case 1, we also consider the model for calibration case 2 to be falsified because the calibrated value of κ=1529.2 d^{−1} is even larger than in case 1 and therefore even less physical.
The compromised fit in calibration case 2 indicated that the effective vertical hydraulic conductivity likely decreases as the thaw depth increases, which would be consistent with observations of reduced hydraulic conductivities at depth. In calibration case 3, we implemented vertical hydraulic conductivity as a function of thaw depth. This model refinement allowed the calibration to achieve a good overall fit to the data (RMSE = 0.37 cm, R^{2}=0.96), as shown in Fig. 3 (red line). K_{r} is well within the estimated limits of 0.7 to 84 m d^{−1} based on a tracer test at a nearby lowcentered polygon conducted by Wales et al. (2020). Based on a rim width of 1 m, the discharge conductance of κ=3.3 d^{−1} would correspond to a drainage boundary layer hydraulic conductivity of 3.3 m d^{−1}, which is also physically reasonable as opposed to the first two calibration estimates.
The standard errors in the calibrated parameters for calibration case 3 listed in Table 1 indicate how well constrained the parameters are by the calibration. It is apparent that the hydraulic conductivities (horizontal and minimum and maximum vertical) are not well constrained with relatively large standard errors. These parameters (or their metaparameters) also have large covariances with each other, indicating their correlated effect on the model. However, despite the lack of constraint of these parameters due to their correlated effect on the model, the calibration does identify reasonable values for them. The standard errors in the discharge conductance, initial polygoncenter water level, and precipitation multiplier indicate that they are well constrained by the calibration.
The calibration verifies that the model is able to capture icewedge polygon drainage characteristics. In the next sections, we perform sensitivity analyses using nondimensional forms of this verified analytical solution to gain insights into icewedge polygon drainage characteristics. The use of nondimensional solution snapshots eliminates the need to consider the precipitation multiplier and thawdepthdependent vertical hydraulic conductivity explicitly. Instead, their effects are implicit in the relative differences between snapshots.
3.2 Drainage pathways
The relative effect of aspect ratio on drainage pathways when the hydraulic conductivities are isotropic (${K}_{\mathrm{r}}/{K}_{\mathrm{z}}=\mathrm{1}$) is shown in Fig. 4. Under isotropic hydraulic conductivities, results show that the drainage pathway for highaspectratio polygons will be predominantly isolated to an annular region at the periphery of the polygon center. As the aspect ratio decreases (in other words, as the thawed region of the polygon subsurface becomes deeper with respect to its width), the portion of the polygon accessed by drainage increases. For an aspect ratio of 20 (Fig. 4a), 95 % of the drainage is focused within around 6 % of the polygon volume, while for an aspect ratio of 2.5 (Fig. 4d), the drainage is spread over around 33 % of the polygon volume. The spreading (increase in the accessed volume) occurs along the radial direction, where the horizontal extent of the accessed volume moves towards the middle of the polygon center (r=0), while the vertical extent is nearly unchanged. The results in Fig. 4 indicate that throughout the thaw season as the thaw depth increases or over successive years as the active layer thickens, the drainage path will spread out towards the middle of the polygon center. Similarly, Fig. 4 can be used to evaluate drainage pathways of polygons of different widths but similar thaw depths. In this context, Fig. 4 indicates that wider polygons will have more focused drainage, while drainage for smaller polygons will be more dispersed.
The effect of anisotropy coefficient on drainage pathways when the aspect ratio is held constant ($R/L=\mathrm{10}$) is shown in Fig. 5. As the anisotropy coefficient increases (i.e., greater preferential horizontal flow), the region accessed by drainage flow becomes larger. Similar to a decreasing aspect ratio (smaller and/or more deeply thawed polygons) in Fig. 4, increasing the anisotropy coefficient leads to a larger radial extent of the accessed region, while the vertical extent is nearly unaffected. When the vertical conductivity is 10 times the horizontal (Fig. 5a), only around 3 % of the polygon volume is accessed by 95 % of the drainage. If horizontal conductivity is 100 times the vertical conductivity (Fig. 5d), around 49 % is accessed. These results indicate that anisotropy has a significant impact on icewedge polygon drainage pathways.
The effect of aspect ratio on drainage pathways when the hydraulic conductivities are highly anisotropic (${K}_{\mathrm{r}}/{K}_{\mathrm{z}}=\mathrm{100}$) is shown in Fig. 6. In this case, contrary to the isotropic case in Fig. 4, as the aspect ratio decreases, the accessed volume generally decreases. However, the largest accessed volume is for an aspect ratio of 10 (Fig. 6b), not 20 (Fig. 6a). This nuance in the dependence of accessed volume on aspect ratio with a high anisotropy coefficient is due to the competing effects of radial extension and vertical contraction of the accessed volume as the aspect ratio decreases. By comparing Figs. 4 and 6, it is also apparent that the volume accessed by drainage is generally larger with a higher anisotropy coefficient.
To gain a global perspective on the trends in drainage pathways with aspect ratio and anisotropy coefficient, Fig. 7 maps the percent of the polygon accessed by 95 % of the drainage as a function of aspect ratio and anisotropy coefficient. This illustrates the overall structure of the combined effect of aspect ratio and anisotropy coefficient on accessed volume (focused versus dispersed drainage). As indicated by the annotated points, the trend in accessed volume with respect to aspect ratio in Fig. 4 is represented along the anisotropy coefficient = 1 transect, while the trend in Fig. 6 is represented by the anisotropy coefficient = 100 transect. The trend in accessed volume with respect to anisotropy coefficient in Fig. 5 along the aspect ratio = 10 transect is likewise indicated. There is a curved, spreading region of high accessed volume with increasing aspect ratio and anisotropy coefficient apparent in the contours of Fig. 7. This curved feature represents the optimal balance of radial and vertical extension of the accessed volume region, as described in the previous paragraph. This structure explains the increasing accessed volume with decreasing aspect ratio in Fig. 4 and the maximum accessed volume at an aspect ratio of 10 in Fig. 6. The drainage flow is most spread out (largest accessed volume) when the aspect ratio is large, and the anisotropy coefficient is high (upper right of Fig. 7). The drainage is the most focused (least accessed volume) when the aspect ratio is large, and the anisotropy coefficient is low (lower right of Fig. 7).
3.3 Ponded water depletion
Depletion curves for the nondimensional ponded water height in the polygon center for various aspect ratios (at fixed high and low anisotropy coefficient) and anisotropies are shown in Fig. 8a, b, and c, respectively. Note that since the depletion of the volume of ponded water is directly related to the ponded water height, the plots in Fig. 8 can be used to obtain either, as indicated by the yaxis labels. As a point of reference between depletion curves, the characteristic time (the time when the height or volume of ponded water reaches $\mathrm{1}/e\approx \mathrm{0.37}$ its initial height or volume) is indicated.
It is apparent from these plots that small, deeply thawed (lower aspect ratio) polygons will drain faster than wide, shallowly thawed (high aspect ratio) polygons, while polygons with a higher anisotropy coefficient will drain faster than those with a lower anisotropy coefficient. It is also apparent that the increase in drainage time with aspect ratio is nearly linear, while for anisotropy coefficient, it is exponential. This is further illustrated in Fig.9, where we plot the trend in characteristic times as a function of aspect ratio for various anisotropies (Fig. 9a) and as a function of anisotropies for various aspect ratios (Fig. 9b). Drainage time (represented by characteristic time) increases in a nearly linear fashion, particularly for high anisotropies, with slight upward curvature increasing for low anisotropies. The drainage time decreases in a nearly log–log fashion (exponentially) with increasing anisotropy coefficient, with an identical trend for different aspect ratios.
A global perspective on drainage timing trends, where characteristic time is mapped as a function of aspect ratio and anisotropy coefficient, is shown in Fig. 10. The fastest (shortest) drainage times are achieved with a low aspect ratio and a high anisotropy coefficient, while the slowest (longest) drainage times are achieved with a high aspect ratio and low anisotropy coefficient. These trends indicate that given the same thawed soil layer thickness, wider polygons will drain more slowly than small polygons. Or, for polygons with similar aspect ratios, those with preferential horizontal flow will drain most quickly, while those with preferential vertical flow will drain most slowly.
3.4 Counteracting effects of aspect ratio and anisotropy
It is important to note that aspect ratio and anisotropy have similar effects on drainage. Within the model used here, in a similar fashion to aspect ratio, the anisotropy coefficient stretches the domain by multiplying the nondimensional radius by $\sqrt{({K}_{\mathrm{z}}/{K}_{\mathrm{r}})}$ (refer to the definition of r^{*} in Eq. 1). As a result, the effect of increasing anisotropy coefficient in the analytical solution is mathematically equivalent to a decrease in aspect ratio. In the end, the equipotential heads and streamlines computed on stretched or compressed radial coordinates (r^{*}) are assigned back to nonmodified radial coordinates (r). This would lead one to believe that it should be possible to obtain two scenarios with different aspect ratios and anisotropies that produce the same mathematical solution. For example, doubling the aspect ratio should produce the mathematically identical effect as dividing the anisotropy coefficient by four. However, this is not the case because the solution involves a Biot parameter (Bi), which defines the ratio of the ability for fluid to conduct across the drainage interface (the vertical outer boundary) relative to the internal hydraulic conductivity as
where κ characterizes the hydraulic conductance across the outer vertical boundary of the model (refer to Appendix A for more details). A large Bi indicates that the outlet boundary interface of the model will not limit the drainage, while a small Bi will result in outletboundaryinterfacelimited drainage. Therefore, to obtain an identical mathematical drainage solution using different combinations of aspect ratio and anisotropy coefficient, one would also need to ensure that Bi is not modified in the process. However, given two solutions, this requires satisfying the conflicting conditions that ${K}_{z\mathrm{1}}/{K}_{r\mathrm{1}}={K}_{z\mathrm{2}}/{K}_{r\mathrm{2}}$ and K_{r1}K_{z1}=K_{r2}K_{z2} (refer to Eqs. 1 and 5), where subscripts 1 and 2 refer to the two solutions. Since these two constraints cannot be simultaneously met, κ or L must also be modified along with aspect ratio and anisotropy coefficient to achieve an identical mathematical solution of drainage. An example where identical mathematical solutions for drainage are obtained is provided in Fig. 11, where the bottom plot is obtained by taking the properties of the top plot and dividing the aspect ratio by 8, dividing the anisotropy coefficient by 8^{2}, and dividing κ by 8. While the solutions are mathematically identical, note that the axes in both plots are not to scale (actual proportions are provided as insets to the plots) and that there are differences in some relative dimensions indicated in the figure (for example, the relative width of the polygon volume accessed by 95 % of the drainage flow is around 4 L in the top case and around 0.5 L in the bottom case).
In the bottom plot in Fig. 11, we demonstrate that despite these differences, the depletion curves are identical for both cases. This result is due to our choice to modify R to change the aspect ratio and modify K_{r} to change the anisotropy coefficient. If either L or K_{z} is chosen to modify the aspect ratio or anisotropy coefficient, respectively, the depletion curves would not necessarily be identical. For example, if L is multiplied by 8, K_{z} multiplied by 8^{2}, and κ divided by 8 (equivalent modifications to aspect ratio and anisotropy coefficient as above), the drainage flow net would still be mathematically equivalent, but the characteristic time would be 8 times shorter than in the original case (refer to Eq. A11). Therefore, while it is possible to obtain mathematically equivalent drainage patterns with counteracting modifications to aspect ratio, anisotropy coefficient, and outflow conductance, the temporal depletion will only be equivalent if R and K_{r} are used to modify the aspect ratio and anisotropy coefficient, respectively.
Our analysis provides new insights into the relative effects of geometry and anisotropy on the manner in which inundated icewedge polygons retain and slowly release water from their centers to their troughs, which form the drainage network of polygonal tundra landscapes. Using a mathematical representation of inundated icewedge polygon drainage (Zlotnik et al., 2020) based on extensive field observations (Helbig et al., 2013; Koch et al., 2014; Liljedahl and Wilson, 2016; Wales et al., 2020) and verified here through calibration to field measurements, we quantify the sensitivity of inundated icewedge polygon drainage to representative polygon sizes, inter and intraannual changes in thaw depth, and preferential flow (hydraulic conductivity anisotropy).
4.1 Calibration implications
The calibration identifies factors which need to be considered by any hydrologic model to simulate drainage from an inundated polygon center. Using a parsimonious model, we were able to identify the refinements required in a polygon drainage model to capture center water levels. Using more complex models would likely obfuscate the identification of these refinements. The final model formulation provides a fast model for predicting the manner and timing of polygon drainage driven by environmental factors. The calibration of the model, driven by transient boundary conditions, provides confidence in the realworld applicability of the sensitivity analysis based on comparisons of steadystate snapshots of the model.
The first refinement is a precipitation multiplier and is based on a simple mass balance indicating that the measured precipitation cannot account for the total increase in ponded water levels after precipitation events. The precipitation multiplier accounts for the fact that precipitation will run off from the rims into the center pond, resulting in a larger increase in ponded water than precipitation measured by a rain gauge. The precipitation multiplier also accounts for any potential rain gauge undercatch. The second refinement is a thawdepthdependent vertical hydraulic conductivity that accounts for the decreased hydraulic conductivity of deeper soil layers. This refinement accounts for the decrease in vertical hydraulic conductivities as the thaw depth increases and includes less hydraulically conductive soils. The final calibration provides a good match to measurements, with a subcentimeter RMSE of 0.37 cm and an R^{2} (equivalent to the Nash–Sutcliffe efficiency here) of 0.96. A key insight from our research is that polygon drainage models need to consider decreases in effective vertical hydraulic conductivity with increasing thaw depth. This research enhances our hydrologic understanding of an Arctic landscape that is undergoing rapid transition due to a warming climate. The hydrology of these landscapes has important implications concerning carbon transport and emissions, subsidence, and Arctic shrubification, to name a few Arctic processes of concern.
Due to covariance in the effects of the hydraulic conductivity parameters on water levels, the hydraulic conductivity parameters are loosely constrained by the calibration (based on local sensitivities). This indicates the importance of constraining the hydraulic conductivity parameters with field measurements if possible (field measurements are not available in our case). However, despite relatively large standard errors in the hydraulic conductivity parameter estimates, the calibration identifies physically realistic values. It should also be noted that despite vertical hydraulic conductivity being loosely constrained in the final calibration (calibration case 3), the model was unable to match water levels with a constant vertical hydraulic conductivity (calibration case 2). This indicates the importance of the thawdepthdependent vertical hydraulic conductivity in the model. The other parameters (discharge conductance, initial polygoncenter water level, and precipitation multiplier) are all well constrained by the analysis.
The water level (Liljedahl and Wilson, 2016) and temperature (Romanovsky et al., 2017) measurements are well constrained with high degrees of resolution and accuracy. While the resolution of precipitation measurements is high (0.1 mm), there is the potential for undercatch during windy precipitation events. More importantly with regard to ponded water levels, the precipitation measurements do not account for the runoff of water during precipitation events from rims into the polygon center, resulting in measured precipitation less than increases in pond water levels. This uncertainty is accounted for through the calibration of a precipitation multiplier, which effectively captures the effect of runoff from the polygon watershed into the center pond along with any potential undercatch during windy precipitation events.
4.2 Analysis implications
A key result of this study is that the geometry and anisotropy of the polygon subsurface have a significant effect on the region of the polygon subsurface predominantly accessed by drainage flow and the time to transition from inundated to drained. Due to the geometry of inundated polygons, the primary drainage pathway is restricted to an annular, radially peripheral region of the icewedge polygon center near the rim. As a result, the middle and lower portions of the polygon center are excluded from the majority of the drainage to varying degrees depending on the polygon geometry and anisotropy. Also, the majority of the ponded water will flow towards the rim of the polygon before infiltrating into the subsurface. Field observations have indicated not only the existence of intrapolygon biogeochemical diversity (Zona et al., 2011; Newman et al., 2015) but also that, based on mineral and nutrient loading, pervasive subsurface flow from centers to troughs exists (Koch et al., 2014). Based on our analysis, polygon geometry and anisotropy will have important effects on the biogeochemistry of polygons (Heikoop et al., 2015; Throckmorton et al., 2015; Newman et al., 2015; Wales et al., 2020; Plaza et al., 2019), affecting dissolved organic matter and mineral flushing from the subsurface of polygon centers to the surface waters of troughs. This is important not only simply concerning the discharge of dissolved organic matter from polygonal tundra landscapes but also for the biogeochemical effects due to photolysis of dissolved organic matter newly exposed to sunlight (Laurion and Mladenov, 2013; Cory et al., 2014).
The potential that the annular, radially peripheral region near the rims will be well flushed of nutrients, while the middle may not, indicates the need for additional field studies designed to measure the effects of anisotropy and preferential flow paths on thermal hydrology and biogeochemistry. For isotropic cases, it should also be considered that the drainage will spread out further towards the middle of the polygon center as the thaw season progresses, and the thawed soil layer thickens (in other words, the aspect ratio decreases; Figs. 4 and 6). However, with a high anisotropy coefficient, the drainage still spreads out towards the polygon middle as the thaw season progresses (aspect ratio decreases), but the depth of the main drainage pathways contracts towards the ground surface (Fig. 6). Ultimately, the effect of this vertical contraction outweighs the lateral spreading, leading to a smaller region being accessed by drainage. Therefore, in the case of smaller and/or more deeply thawed (low aspect ratio) polygons with a high anisotropy coefficient, the region close to the surface will be flushed more than the deeper regions. Deeper nutrients that become available for transport because of recent thaw may still be less accessible to aqueous transport due to drainage pathways in polygons with a high anisotropy coefficient. Field measurements of geochemical depth profiles in polygons support limited deep flushing. As reviewed in Newman et al. (2015), concentration increases in geogenic solutes and organic carbon with activelayer depth are relatively common in polygonal terrain. It is unlikely that these sometimes large geochemicaldepth gradients would be preserved with significant amounts of deep flushing. Interannual deepening of the active layer will lead to even more dramatic evolution in drainage patterns described above during the thaw season as the range of thaw depths encountered increases. This research reinforces the need for field studies on anisotropy and preferential flow in polygon landscapes to better understand the hydrologic transitions and feedbacks that will occur in a warming climate.
For a given thaw depth, advective heat transport will be more focused near the rim for larger polygons and may result in enhanced icewedge degradation (Wright et al., 2009). Based on the relative drainage pathways, advective heat transport will be most pronounced in larger and/or shallower (large aspect ratio) polygons with a low anisotropy coefficient (refer to Figs. 4a and 5a). Therefore, drainage events for wide, isotropic (or with preferential vertical flow) polygons may result in enhanced icewedge top thawing, which would promote low to highcentered polygon transition. The more spread out the drainage is throughout the thawed soil layer, the less pronounced this effect may be (for example, wide polygons with a high anisotropy coefficient as in Fig. 6a and b). However, small polygons with a high anisotropy coefficient will restrict the drainage flow near the ground surface, potentially reducing the advective transport of heat to the icewedge top as well (Fig. 6d). Similarly, increasing seasonal thaw depth and interannual activelayer thickness will lead to polygons with less focused advective heat transport towards icewedge tops, potentially providing a negative feedback on icewedge degradation. These results provide an additional perspective on icewedge degradation, complementing previous research that found that the process is also strongly controlled by geometry (Abolt et al., 2020) and hydrologic conditions (Nitzbon et al., 2019).
Small polygons with deeply thawed soil layers (low aspect ratios) and high horizontal preferential flow (high anisotropy coefficient) have the potential to drain most quickly. Therefore, all other factors being equal, regions of polygonal tundra characterized by small, deeply thawed, anisotropic polygons will drain more quickly and consequently will have a greater potential for nutrient flushing, transition from methane to carbon dioxide atmospheric emissions, and biological succession than regions with large, shallowly thawed, isotropic polygons. Concerning temporal changes during the thaw season, polygon pond depletion will slow down as the thaw depth increases, and this reduction will become more dramatic over the thaw season as interannual activelayer thickness increases. This implies that a thickening active layer due to warming trends may result in slower pond depletion. For a given location, factors such as regional flow patterns, largescale topography, etc. will influence the region's overall drainage timing. However, along with these other factors, our analysis indicates that aspect ratio will have a nearly linear positive relationship, while the anisotropy coefficient will have an exponential negative relationship with drainage timing.
4.3 Model limitations
The drainage pathways and timing presented here are based on hydrogeological first principles (Harr, 1962; Cedergren, 1968; Freeze and Cherry, 1979; Bear, 1979) using a generalization of icewedge polygon geometry and hydraulic properties that allow insights into drainage pathways and timing. The model captures the physical forces that ponded water exerts on a cylindrical porous disc with drainage allowed radially through its sides. While the cylindrical geometry and anisotropy coefficients considered here do not cover all potential variations present in icewedge polygons, the impact of those variations will cause deviations from the idealized scenarios considered here. For example, heterogeneities will warp the base case hydraulic head equipotentials and streamlines. The cylindrical idealization of icewedge polygon geometry will best approximate drainage in nearly symmetrical hexagonal icewedge polygons, while nonsymmetric or square icewedge polygon drainage will deviate most from that shown here. The model assumes uniform radial discharge around the periphery of the polygon center. Low points in the elevated frozen ground under the rim will lead to deviations from our results. Elevated frozen ground under the rim will also warp the streamlines upwards within a localized region near the rim, but this effect will dissipate a short distance into the polygon center.
In our sensitivity analysis, we use idealized conceptualizations and dimensionless variables, which allow hydrologic characteristics of drainage to be exposed in their most fundamental form. To clearly and concisely expose these characteristics, we neglect factors such as evaporation, precipitation, and nonidealized polygon geometry (evaporation and precipitation are included in the calibration). As a result, our analysis is not intended to provide predictive capability across all polygonal tundra scenarios but to provide hydrologic intuition into the relative effects of geometry and anisotropy on inundated icewedge polygon drainage and timing.
As the analytical solution applies to ponded conditions, it does not apply to freezeup at the end of the thaw season, after the ponded water in the polygon center freezes. At this time, the active layer freezes simultaneously from the top and bottom, and cryosuction draws water towards the freezing fronts. This redistribution of water affects icewedge polygon drainage, and Wales et al. (2020) postulate that it may explain some of the observations of their tracer test. More complex models than applied here are required to capture the details of water redistribution during freezeup (Painter, 2011; Atchley et al., 2016; Schuh et al., 2017).
Since the model is based on the saturated groundwater flow equation, in its current form it cannot be applied to noninundated lowcentered polygons. Since it is based on having a ponded center, the model is also not applicable to highcentered polygons. However, a similar approach to that presented here could be taken with the unsaturated groundwater flow equation to capture these other polygon scenarios and types.
4.4 Relation to other mathematical analyses
As our analysis is the first detailed examination of inundated icewedge polygon drainage patterns, our results provide a new perspective to existing mathematical analyses. Much of the existing mathematical analyses investigate drainage from polygonal landscapes. For example, CrestoAleina et al. (2013) analyze lowcentered polygonal landscape hydrology stochastically using percolation theory, a method that is not designed to investigate drainage pathways of individual polygons. Nitzbon et al. (2019, 2020) investigate the hydrology of polygonal tundra using a “tiled” approach that is also not designed to map individual icewedge polygon drainage pathways. In these cases, our analysis enhances our understanding by considering details that are glossed over in these analyses for the sake of considering drainage from many polygons. Other research has focused on processrich 1D analyses that are unable to consider lateral flow. For example, Atchley et al. (2015), Harp et al. (2016), and Atchley et al. (2016) use a processrich, complex 1D model to calibrate and perform uncertainty and sensitivity analyses of thermalhydrology models of icewedge polygons. Our analysis augments these by providing insights into lateral fluid flow. Abolt et al. (2020) perform a sensitivity analysis of the thermal effects of icewedge thaw on individual polygons using a fullphysics thermalhydrology model. Our analysis adds to this research by providing indications of the drainage patterns that would exist along with the thermal effects. Jan et al. (2020) calibrate a fullphysics thermalhydrology model of icewedge polygons using field data, but since the focus of their research is calibration, they do not investigate drainage pathways. As a result, our analysis provides new information that augments existing research, helps the Arctic hydrology research community gain an intuitive understanding of inundated icewedge polygon drainage, and provides new modeling and experimental research directions based on the idealized base cases considered here.
Our results affirm that existing conceptualizations of polygon drainage need revisiting, as was implied by a preliminary analysis conducted by Zlotnik et al. (2020). Fundamental hydrological principles indicate that, due to lowcentered polygon geometry and hydraulic properties, inundated polygon drainage flux will not be uniform throughout the thawed soil layer. This result indicates that while some portions of the polygon thaw layer will be well flushed, other portions will be poorly flushed. This has important implications for the aqueous transport of carbon and other dissolved nutrients from polygonal tundra. It also indicates the potential for focused advective heat transport from the ponded water towards the icewedge tops, which could enhance icewedge degradation. Within this context, we provide the following conclusions:

A simple analytical solution based on hydrological first principles is able to capture icewedge polygon drainage dynamics over an entire thaw season. Consideration within the model of runoff from topographic highs (rims) during rain events and decreasing effective vertical hydraulic conductivity with increasing thaw depth is required to match field measurements. We were able to identify these necessary model components by using a parsimonious model that has the ability to fail.

We provide rigorous confirmation that the majority of drainage from inundated icewedge polygon centers occurs along an annular region along their radial periphery; however, polygon geometry and hydraulic conductivity anisotropy significantly impact the drainage pathways, as originally postulated by Zlotnik et al. (2020). This implies that nutrient flushing and the advective thermal transfer will be focused along polygon edges, neglecting the polygon center. Additionally, our results indicate the following:

A combination of high aspect ratio (wide, shallow polygons) and high anisotropy coefficient (preferential horizontal flow) results in the greatest spreading of drainage flow towards the middle of the polygon center and the largest fraction of the polygon volume being accessed by drainage flow. In these cases, the nutrient flushing will be more uniform than other cases, and advective heat transport towards the icewedge top will be less focused and therefore less able to thaw the icewedge top.

A combination of high aspect ratio (wide, shallow polygons) and low anisotropy coefficient (preferential vertical flow) results in the greatest focusing of drainage flow and the smallest fraction of the polygon volume being accessed by drainage flow. In these cases, nutrient flushing will be more localized along the outer edge of the polygon (nonuniform) and the advective heat transport towards the icewedge tops most focused, possibly resulting in greater degradation of ice wedges than in other cases.

Combinations of aspect ratio and anisotropy coefficient have counteracting effects of radial versus vertical extension and contraction of drainage pathways, producing nonmonotonic relationships between aspect ratio or anisotropy coefficient and accessed volume (curved feature in accessed volume response surface in Fig. 7). Therefore, the combined nonlinear effects of geometry and anisotropy on drainage patterns must be considered when evaluating nutrient flushing and advective heat transport.


Polygon drainage time has an approximately positive linear relationship with aspect ratio when anisotropy is held constant; in other words, wide, shallow polygons drain slowly, while small, deep polygons drain quickly. This implies that polygonal tundra with larger polygons may drain more slowly than tundra composed of smaller polygons.

Polygon drainage time has a negative exponential relationship with anisotropy coefficient when aspect ratio is held constant. In other words, increases in preferential horizontal flow lead to exponentially faster drainage. Therefore, polygonal tundra with greater preferential horizontal flow, due to more pronounced horizontal stratigraphy or ice lens development, will drain faster, while less horizontally stratified tundra, due to cryoturbation, for example, will drain more slowly.
Here, we present the analytical solutions for hydraulic heads and stream function under the center of an inundated icewedge polygon and the depletion curve of the ponded water height due to drainage to the surrounding trough. For details on the derivations and validation of the model, refer to Zlotnik et al. (2020).
A1 Hydraulic head and stream function analytical solutions
We idealize the subsurface below the center region of a lowcentered polygon as a thin cylinder with its radius in the horizontal direction and length in the vertical direction (refer to Fig. 1). Based on this approximation, the hydraulic heads ($h(r,z,t)$) in an inundated polygon will satisfy the following equation:
However, considering that storativity is negligible given the vertical and horizontal scale of the domain, the righthand term can be neglected as
where h is the hydraulic head, r is the radial coordinate, z is the depth coordinate (positive in the downward vertical direction), K_{r} is the horizontal (radial) hydraulic conductivity, K_{z} is the vertical hydraulic conductivity, t is time, and S_{s} is specific storage. The system is fully specified by ensuring that (BC1) the heads along the central vertical axis ($h(\mathrm{0},z,t)$) are finite, (BC2) the heads along the ground surface are equal to the height of ponded water in the polygon center ($h(r,\mathrm{0},t)={H}_{\mathrm{c}}\left(t\right)$), (BC3) the change in heads along the outer vertical boundary is governed by the change in heads along the boundary and the height of water in the trough (H_{t})
(BC4) the bottom of the model has a zero head gradient
and (BC5) the change in volume of ponded water in the polygon center is related to the vertical head gradients along the ground surface
where R is the radius of the polygon, L is the depth of the thawed subsurface of the polygon, H_{c}(t) is the ponded water height in the center of the polygon at time t, H_{t} is the height of water in the trough, and κ characterizes the conductance to flow across the drainage interface to the trough (in our case, the hydraulic conductance of the soil layer under the rim).
Using dimensionless coordinates and parameters defined as
dimensionless solutions for hydraulic heads and the stream function can be obtained as
and
respectively, where J_{m}, with m=0 or 1, is the Bessel function of the first kind of mth order, and λ_{n} is the nth root of the equation
The solutions can be verified by direct substitution of Eqs. (A7) and (A8) into the boundary value problem defined by Eqs. (A2) to (A6).
A2 Ponded height depletion curve analytical solution
The ponded height depletion curve can be defined as
where t_{L} is the characteristic depletion time defined as
which is the time when ${H}_{\mathrm{c}}\left(t\right)/{H}_{\mathrm{c},\mathrm{0}}=\mathrm{1}/e$, in other words, when the ponded height is approximately 37 % of its initial height. The function $F(\mathrm{Bi},{R}^{*})$ can be evaluated as
The depletion curve can be expressed in nondimensional terms as
where nondimensional time ${t}^{*}=t/{t}_{L}$. The solution can be verified by direct substitution of Eqs. (A10) into (A5).
MATLAB code of the analytical solution is available via Zlotnik et al. (2020) at http://www.mdpi.com/20734441/12/12/3376/s1, last access: 18 August 2021.
All data sets are freely available to the public at the provided references in the paper.
A video illustrating the validation of the analytical solution to an inundated icewedge polygon drainage event in 2012 near Utqiaġvik is available via Zlotnik et al. (2020) at http://www.mdpi.com/20734441/12/12/3376/s1, last access: 18 August 2021.
DRH and VZ developed the conceptual model of inundated polygonal tundra hydrology. VZ derived the analytical solutions. VZ and DRH encoded the analytical solutions. CJA created Fig. 1. DRH created Figs. 4–11. ALA provided text for the “Introduction” and “Discussion” sections. EJ provided text for the discussion section. ALA, BDN, EJ, and CJW provided critical reviews of the manuscript. CJW secured funding.
The authors declare that they have no conflict of interest.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
The Next Generation Ecosystem Experiments Arctic (NGEEArctic) project (DOE ERKP757), funded by the Office of Biological and Environmental Research within the US Department of Energy's Office of Science, supported this research. Sofia Avendaño and Bulbul Ahmmed provided reviews during the development of this article, providing technical and editorial improvements.
This research has been supported by the US Department of Energy, Office of Science (grant no. DOE ERKP757).
This paper was edited by Moritz Langer and reviewed by Jan Nitzbon and one anonymous referee.
Abnizova, A., Siemens, J., Langer, M., and Boike, J.: Small ponds with major impact: The relevance of ponds and lakes in permafrost landscapes to carbon dioxide emissions, Global Biogeochem. Cy., 26, 1–9, 2012. a
Abolt, C. J. and Young. M.H.: Highresolution mapping of spatial heterogeneity in ice wedge polygon geomorphology near Prudhoe Bay, Alaska, Scientific data, 7, 1–7, 2020. a, b
Abolt, C. J., Young, M. H., Atchley, A. L., and Harp, D. R.: Microtopographic control on the ground thermal regime in ice wedge polygons, The Cryosphere, 12, 1957–1968, https://doi.org/10.5194/tc1219572018, 2018. a
Abolt, C. J., Young, M. H., Atchley, A. L., and Wilson, C. J.: Brief communication: Rapid machinelearningbased extraction and measurement of ice wedge polygons in highresolution digital elevation models, The Cryosphere, 13, 237–245, https://doi.org/10.5194/tc132372019, 2019. a
Abolt, C. J., Young, M. H., Atchley, A. L., Harp, D. R., and Coon, E. T.: Feedbacks between surface deformation and permafrost degradation in ice wedge polygons, Arctic Coastal Plain, Alaska, J. Geophys. Res.Earth, 125, e2019JF005349, https://doi.org/10.1029/2019JF005349, 2020. a, b
Atchley, A. L., Painter, S. L., Harp, D. R., Coon, E. T., Wilson, C. J., Liljedahl, A. K., and Romanovsky, V. E.: Using field observations to inform thermal hydrology models of permafrost dynamics with ATS (v0.83), Geosci. Model Dev., 8, 2701–2722, https://doi.org/10.5194/gmd827012015, 2015. a, b
Atchley, A. L., Coon, E. T., Painter, S. L., Harp, D. R., and Wilson, C. J.: Influences and interactions of inundation, peat, and snow on active layer thickness, Geophys. Res. Lett., 43, 5116–5123, 2016. a, b, c
Bear, J.: Hydraulics of groundwater New York, Mc GrawHill Inc, New York, 1979. a, b, c
Beckwith, C. W., Baird, A. J., and Heathwaite, A. L.: Anisotropy and depthrelated heterogeneity of hydraulic conductivity in a bog peat. I: laboratory measurements, Hydrol. Process., 17, 89–101, 2003. a, b
Billings, W. and Peterson, K.: Vegetational change and icewedge polygons through the thawlake cycle in Arctic Alaska, Arct. Alp. Res., 12, 413–432, 1980. a
Boike, J., Wille, C., and Abnizova, A.: Climatology and summer energy and water balance of polygonal tundra in the Lena River Delta, Siberia, J. Geophys. Res.Biogeo., 113, 1–15, 2008. a, b, c
Brown, J., Ferrians Jr, O., Heginbottom, J., and Melnikov, E.: CircumArctic map of permafrost and groundice conditions, US Geological Survey Reston, VA, 1997. a
Brown, V. A., McDonnell, J. J., Burns, D. A., and Kendall, C.: The role of event water, a rapid shallow flow component, and catchment size in summer stormflow, J. Hydrol., 217, 171–190, 1999. a
Cedergren, H. R.: Seepage, drainage, and flow nets, John Wiley & Sons, New York, 1968. a, b
Conway, T. and Steele, L.: Carbon dioxide and methane in the Arctic atmosphere, J. Atmos. Chem., 9, 81–99, 1989. a
Cory, R. M., Ward, C. P., Crump, B. C., and Kling, G. W.: Sunlight controls water column processing of carbon in arctic fresh waters, Science, 345, 925–928, 2014. a
Cresto Aleina, F., Brovkin, V., Muster, S., Boike, J., Kutzbach, L., Sachs, T., and Zuyev, S.: A stochastic model for the polygonal tundra based on Poisson–Voronoi diagrams, Earth Syst. Dynam., 4, 187–198, https://doi.org/10.5194/esd41872013, 2013. a, b
Freeze, R. A. and Cherry, J. A.: Groundwater, PrenticeHall, Englewood Cliffs, New Jersey, 1979. a, b
Harp, D. R., Atchley, A. L., Painter, S. L., Coon, E. T., Wilson, C. J., Romanovsky, V. E., and Rowland, J. C.: Effect of soil property uncertainties on permafrost thaw projections: a calibrationconstrained analysis, The Cryosphere, 10, 341–358, https://doi.org/10.5194/tc103412016, 2016. a, b
Harr, M.: Groundwater and Seepage, McGrawHill, New York, 1962. a, b
Heikoop, J. M., Throckmorton, H. M., Newman, B. D., Perkins, G. B., Iversen, C. M., Roy Chowdhury, T., Romanovsky, V., Graham, D. E., Norby, R. J., Wilson, C. J., and Wullschleger, S. D.: Isotopic identification of soil and permafrost nitrate sources in an Arctic tundra ecosystem, J. Geophys. Res.Biogeo., 120, 1000–1017, 2015. a, b
Helbig, M., Boike, J., Langer, M., Schreiber, P., Runkle, B. R., and Kutzbach, L.: Spatial and seasonal variability of polygonal tundra water balance: Lena River Delta, northern Siberia (Russia), Hydrogeol. J., 21, 133–147, 2013. a, b, c, d, e, f, g, h, i
Hinzman, L., Busey, B., Cable, W., and Romanovsky, V.: Surface meteorology, Barrow, Alaska, Area A, B, C and D, ongoing from 2012, Next Generation Ecosystem Experiments Arctic Data Collection, Oak Ridge National Laboratory, U.S. Department of Energy, Oak Ridge, Tennessee, USA, Data accessed on 28 December 2020, https://doi.org/10.5440/1164893, 2014. a
Hugelius, G., Strauss, J., Zubrzycki, S., Harden, J. W., Schuur, E. A. G., Ping, C.L., Schirrmeister, L., Grosse, G., Michaelson, G. J., Koven, C. D., O'Donnell, J. A., Elberling, B., Mishra, U., Camill, P., Yu, Z., Palmtag, J., and Kuhry, P.: Estimated stocks of circumpolar permafrost carbon with quantified uncertainty ranges and identified data gaps, Biogeosciences, 11, 6573–6593, https://doi.org/10.5194/bg1165732014, 2014. a
Jafarov, E., Parsekian, A., Schaefer, K., Liu, L., Chen, A., Panda, S., and Zhang, T.: Estimating active layer thickness and volumetric water content from ground penetrating radar measurements in Barrow, Alaska, Geosci. Data J., 4, 72–79, 2017. a
Jan, A., Coon, E. T., Graham, J. D., and Painter, S. L.: A subgrid approach for modeling microtopography effects on overland flow, Water Resour. Res., 54, 6153–6167, 2018a. a
Jan, A., Coon, E. T., Painter, S. L., Garimella, R., and Moulton, J. D.: An intermediatescale model for thermal hydrology in lowrelief permafrostaffected landscapes, Computat. Geosci., 22, 163–177, 2018b. a
Jan, A., Coon, E. T., and Painter, S. L.: Evaluating integrated surface/subsurface permafrost thermal hydrology models in ATS (v0.88) against observations from a polygonal tundra site, Geosci. Model Dev., 13, 2259–2276, 2020. a
Jorgenson, M. T., Shur, Y. L., and Pullman, E. R.: Abrupt increase in permafrost degradation in Arctic Alaska, Geophys. Res. Lett., 33, 1–4, 2006. a
Jorgenson, M. T., Kanevskiy, M., Shur, Y., Moskalenko, N., Brown, D., Wickland, K., Striegl, R., and Koch, J.: Role of ground ice dynamics and ecological feedbacks in recent ice wedge degradation and stabilization, J. Geophys. Res.Earth, 120, 2280–2297, 2015. a
King, T. V., Neilson, B. T., Overbeck, L. D., and Kane, D. L.: A distributed analysis of lateral inflows in an Alaskan Arctic watershed underlain by continuous permafrost, Hydrol. Process., 34, 633–648, 2020. a
Koch, J. C.: Lateral and subsurface flows impact arctic coastal plain lake water budgets, Hydrol. Process., 30, 3918–3931, 2016. a, b, c, d
Koch, J. C., Gurney, K., and Wipfli, M. S.: Morphologydependent water budgets and nutrient fluxes in Arctic thaw ponds, Permafrost Periglac., 25, 79–93, 2014. a, b, c, d, e, f, g
Koch, J. C., Jorgenson, M. T., Wickland, K. P., Kanevskiy, M., and Striegl, R.: Ice wedge degradation and stabilization impact water budgets and nutrient cycling in Arctic trough ponds, J. Geophys. Res.Biogeo., 123, 2604–2616, 2018. a, b, c, d, e
Lachenbruch, A. H.: Mechanics of thermal contraction cracks and icewedge polygons in permafrost, Vol. 70, Geol. Soc. Am., 70, 1–63, 1962. a
Lara, M. J., McGuire, A. D., Euskirchen, E. S., Tweedie, C. E., Hinkel, K. M., Skurikhin, A. N., Romanovsky, V. E., Grosse, G., Bolton, W. R., and Genet, H.: Polygonal tundra geomorphological change in response to warming alters future CO_{2} and CH_{4} flux on the Barrow Peninsula, Glob. Change Biol., 21, 1634–1651, 2015. a
Larouche, J. R., Abbott, B. W., Bowden, W. B., and Jones, J. B.: The role of watershed characteristics, permafrost thaw, and wildfire on dissolved organic carbon biodegradability and water chemistry in Arctic headwater streams, Biogeosciences, 12, 4221–4233, https://doi.org/10.5194/bg1242212015, 2015. a
Laurion, I. and Mladenov, N.: Dissolved organic matter photolysis in Canadian arctic thaw ponds, Environ. Res. Lett., 8, 035026, https://doi.org/10.1088/17489326/8/3/035026, 2013. a, b
Leffingwell, E. d. K.: Groundice wedges: The dominant form of groundice on the north coast of Alaska, J. Geol., 23, 635–654, 1915. a, b
Lewkowicz, A. G.: Icewedge rejuvenation, fosheim peninsula, ellesmere Island, Canada, Permafrost Periglac., 5, 251–268, 1994. a
Liljedahl, A. K. and Wilson, C. J.: Ground Water Levels for NGEE Areas A, B, C and D, Barrow, Alaska, 2012–2014, Next Generation Ecosystem Experiments Arctic Data Collection, Oak Ridge National Laboratory, U.S. Department of Energy, Oak Ridge, Tennessee, USA, Data accessed on 8 February 2020, https://doi.org/10.5440/1183767, 2016. a, b, c, d
Liljedahl, A. K., Boike, J., Daanen, R. P., Fedorov, A. N., Frost, G. V., Grosse, G., Hinzman, L. D., Iijma, Y., Jorgenson, J. C., Matveyeva, N., and Necsoiu, M.: PanArctic icewedge degradation in warming permafrost and its influence on tundra hydrology, Nat. Geosci., 9, 312–318, 2016. a, b, c, d
Mackay, J.: Thermally induced movements in icewedge polygons, western Arctic coast: a longterm study, Geogr. Phys. Quatern., 54, 41–68, 2000. a, b, c
Mackay, J. R.: The world of underground ice, Ann. Assoc. Am. Geogr., 62, 1–22, 1972. a
Mackay, J. R.: Active layer slope movement in a continuous permafrost environment, Garry Island, Northwest Territories, Canada, Can. J. Earth Sci., 18, 1666–1680, 1981. a
Mackay, J. R.: Some observations on the growth and deformation of epigenetic, syngenetic and antisyngenetic ice wedges, Permafrost Periglac., 1, 15–29, 1990. a
Matsuoka, N. and Moriwaki, K.: Frost heave and creep in the Sør Rondane Mountains, Antarct. Arct. Alp. Res., 24, 271–280, 1992. a
McDonnell, J., Owens, I. F., and Stewart, M.: A case study of shallow flow paths in a steep zeroorder basin 1, J. Am. Water Resour. Assoc., 27, 679–685, 1991. a
Minke, M., Donner, N., Karpov, N. S., de Klerk, P., and Joosten, H.: Distribution, diversity, development and dynamics of polygon mires: examples from Northeast Yakutia (Siberia), Peatlands Internation, 36–40, 2007. a
Moore, T. and Dalva, M.: Methane and carbon dioxide exchange potentials of peat soils in aerobic and anaerobic laboratory incubations, Soil Biol. Biochem., 29, 1157–1164, 1997. a
Newman, B. D., Throckmorton, H. M., Graham, D. E., Gu, B., Hubbard, S. S., Liang, L., Wu, Y., Heikoop, J. M., Herndon, E. M., Phelps, T. J., Wilson, C. J., and Wullschleger, S. D.: Microtopographic and depth controls on active layer chemistry in Arctic polygonal ground, Geophys. Res. Lett., 42, 1808–1817, 2015. a, b, c, d, e
Nitzbon, J., Langer, M., Westermann, S., Martin, L., Aas, K. S., and Boike, J.: Pathways of icewedge degradation in polygonal tundra under different hydrological conditions, The Cryosphere, 13, 1089–1123, https://doi.org/10.5194/tc1310892019, 2019. a, b, c, d
Nitzbon, J., Westermann, S., Langer, M., Martin, L. C., Strauss, J., Laboor, S., and Boike, J.: Fast response of cold icerich permafrost in northeast Siberia to a warming climate, Nat. Commun., 11, 1–11, 2020. a, b
Norby, R. J., Sloan, V. L., Iversen, C. M., and Childs, J.: Controls on finescale spatial and temporal variability of plantavailable inorganic nitrogen in a polygonal tundra landscape, Ecosystems, 22, 528–543, 2019. a
O'Shea, S. J., Allen, G., Gallagher, M. W., Bower, K., Illingworth, S. M., Muller, J. B. A., Jones, B. T., Percival, C. J., Bauguitte, S. J.B., Cain, M., Warwick, N., Quiquet, A., Skiba, U., Drewer, J., Dinsmore, K., Nisbet, E. G., Lowry, D., Fisher, R. E., France, J. L., Aurela, M., Lohila, A., Hayman, G., George, C., Clark, D. B., Manning, A. J., Friend, A. D., and Pyle, J.: Methane and carbon dioxide fluxes and their regional scalability for the European Arctic wetlands during the MAMM project in summer 2012, Atmos. Chem. Phys., 14, 13159–13174, https://doi.org/10.5194/acp14131592014, 2014. a
Painter, S. L.: Threephase numerical model of water migration in partially frozen geological media: model formulation, validation, and applications, Comput. Geosci., 15, 69–85, 2011. a
Plaza, C., Pegoraro, E., Bracho, R., Celis, G., Crummer, K. G., Hutchings, J. A., Pries, C. E. H., Mauritz, M., Natali, S. M., Salmon, V. G., and Schädel, C.: Direct observation of permafrost degradation and rapid soil carbon loss in tundra, Nat. Geosci., 12, 627–631, 2019. a, b
Pollock, M. D., O'Donnell, G., Quinn, P., Dutton, M., Black, A., Wilkinson, M. E., Colli, M., Stagnaro, M., Lanza, L. G., Lewis, E., and Kilsby, C. G.: Quantifying and mitigating windinduced undercatch in rainfall measurements, Water Resour. Res., 54, 3863–3875, 2018. a
Popper, K.: The logic of scientific discovery, Routledge, 2005. a
Quinton, W. and Marsh, P.: A conceptual framework for runoff generation in a permafrost environment, Hydrol. Process., 13, 2563–2581, 1999. a
Raymond, P. A., Hartmann, J., Lauerwald, R., Sobek, S., McDonald, C., Hoover, M., Butman, D., Striegl, R., Mayorga, E., Humborg, C., and Kortelainen, P.: Global carbon dioxide emissions from inland waters, Nature, 503, 355–359, 2013. a
Raynolds, M. K., Walker, D. A., Ambrosius, K. J., Brown, J., Everett, K. R., Kanevskiy, M., Kofinas, G. P., Romanovsky, V. E., Shur, Y., and Webber, P. J.: Cumulative geoecological effects of 62 years of infrastructure and climate change in icerich permafrost landscapes, Prudhoe Bay Oilfield, Alaska, Glob. Change Biol., 20, 1211–1224, 2014. a
Rodell, M., Houser, P. R., Jambor, U. E. A., Gottschalck, J., Mitchell, K., Meng, C. J., Arsenault, K., Cosgrove, B., Radakovich, J., Bosilovich, M., and Entin, J. K.: The global land data assimilation system, Bull. Am. Meteorol. Soc., 85, 381–394, 2004. a
Romanovsky, V., Cable, W., and Dolgikh, K.: Subsurface temperature, moisture, thermal conductivity and heat flux, Barrow, Area A, B, C, D, Next Generation Ecosystem Experiments Arctic Data Collection, Oak Ridge National Laboratory, U.S. Department of Energy, Oak Ridge, Tennessee, USA, Data accessed on 29 December 2020, https://doi.org/10.5440/1126515, 2017. a, b
Schuh, C., Frampton, A., and Christiansen, H. H.: Soil moisture redistribution and its effect on interannual active layer temperature and thickness variations in a dry loess terrace in Adventdalen, Svalbard, The Cryosphere, 11, 635–651, https://doi.org/10.5194/tc116352017, 2017. a
Schuur, E. A., Bockheim, J., Canadell, J. G., Euskirchen, E., Field, C. B., Goryachkin, S. V., Hagemann, S., Kuhry, P., Lafleur, P. M., Lee, H., and Mazhitova, G.: Vulnerability of permafrost carbon to climate change: Implications for the global carbon cycle, BioScience, 58, 701–714, 2008. a
Shiklomanov, N. I., Streletskiy, D. A., Nelson, F. E., Hollister, R. D., Romanovsky, V. E., Tweedie, C. E., Bockheim, J. G., and Brown, J.: Decadal variations of activelayer thickness in moisturecontrolled landscapes, Barrow, Alaska, J. Geophys. Res.Biogeo., 115, 1–14, 2010. a
Tarnocai, C., Canadell, J., Schuur, E. A., Kuhry, P., Mazhitova, G., and Zimov, S.: Soil organic carbon pools in the northern circumpolar permafrost region, Global Biogeochem. Cy., 23, 1–11, 2009. a
Throckmorton, H. M., Heikoop, J. M., Newman, B. D., Altmann, G. L., Conrad, M. S., Muss, J. D., Perkins, G. B., Smith, L. J., Torn, M. S., Wullschleger, S. D., and Wilson, C. J.: Pathways and transformations of dissolved methane and dissolved inorganic carbon in Arctic tundra watersheds: Evidence from analysis of stable isotopes, Global Biogeochem. Cy., 29, 1893–1910, 2015. a, b
Transtrum, M. K., Machta, B. B., and Sethna, J. P.: Geometry of nonlinear least squares with applications to sloppy models and optimization, Phys. Rev. E, 83, 036701, https://doi.org/10.1103/PhysRevE.83.036701, 2011. a
Vaughn, L. J., Conrad, M. E., Bill, M., and Torn, M. S.: Isotopic insights into methane production, oxidation, and emissions in Arctic polygon tundra, Glob. Change Biol., 22, 3487–3502, 2016. a
Wainwright, H. M., Dafflon, B., Smith, L. J., Hahn, M. S., Curtis, J. B., Wu, Y., Ulrich, C., Peterson, J. E., Torn, M. S., and Hubbard, S. S.: Identifying multiscale zonation and assessing the relative importance of polygon geomorphology on carbon fluxes in an Arctic tundra ecosystem, J. Geophys. Res.Biogeo., 120, 788–808, 2015. a, b
Wales, N. A., GomezVelez, J. D., Newman, B. D., Wilson, C. J., Dafflon, B., Kneafsey, T. J., Soom, F., and Wullschleger, S. D.: Understanding the relative importance of vertical and horizontal flow in icewedge polygons, Hydrol. Earth Syst. Sci., 24, 1109–1129, https://doi.org/10.5194/hess2411092020, 2020. a, b, c, d, e, f, g, h, i, j, k, l
Wolter, J., Lantuit, H., Fritz, M., MaciasFauria, M., MyersSmith, I., and Herzschuh, U.: Vegetation composition and shrub extent on the Yukon coast, Canada, are strongly linked to icewedge polygon degradation, Polar Res., 35, 27489, https://doi.org/10.3402/polar.v35.27489, 2016. a
Wright, N., Hayashi, M., and Quinton, W. L.: Spatial and temporal variations in active layer thawing and their implication on runoff generation in peatcovered permafrost terrain, Water Resour. Res., 45, 1–13, 2009. a
Zhu, X., Zhuang, Q., Gao, X., Sokolov, A., and Schlosser, C. A.: PanArctic land–atmospheric fluxes of methane and carbon dioxide in response to climate change over the 21st century, Environ. Res. Lett., 8, 045003, https://doi.org/10.1088/17489326/8/4/045003, 2013. a
Zlotnik, V. A., Harp, D. R., Jafarov, E. E., and Abolt, C. J.: A Model of Ice Wedge Polygon Drainage in Changing Arctic Terrain, Water, 12, 3376, https://doi.org/10.3390/w12123376, 2020. a, b, c, d, e, f, g, h, i, j, k, l
Zona, D., Lipson, D., Zulueta, R., Oberbauer, S., and Oechel, W.: Microtopographic controls on ecosystem functioning in the Arctic Coastal Plain, J. Geophys. Res.Biogeo., 116, 1–12, 2011. a, b, c, d, e