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 low-centered 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 non-dimensional form of the solution for hydraulic heads and stream function (refer to Eqs. A7 and A8). The non-dimensional 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 non-dimensional 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.

Non-dimensional 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 non-dimensional 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 non-dimensional 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 third-type 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 non-dimensional stream function >0.05).

The change in non-dimensional 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 least-squares calibrations using a Levenberg–Marquardt approach (Transtrum et al., 2011). The objective function is the sum-of-squared errors (SSEs) between the measured and simulated polygon-center water levels expressed as

where θ is a vector of model parameters, and $H_{\mathrm{c}}^m\left(t_i\right)$ is the measured polygon-center 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 depth-dependent 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 meta-parameter ${\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}+\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 super-elliptical function of thaw depth to calibration case 2 as $K_{\mathrm{z}}\left(D\right)=K_{\mathrm{z}}^{\min }+(K_{\mathrm{z}}^{\max }-K_{\mathrm{z}}^{\min }){\left(\mathrm{1}-D_{\mathrm{norm}}^a\right)}^{\mathrm{1}/a},D_{\mathrm{norm}}=\frac{D-D^{\min }}{D^{\max }-D^{\min }}$, where $K_{\mathrm{z}}^{\max }$ and $K_{\mathrm{z}}^{\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 super-elliptical shape parameter that controls curvature. $K_{\mathrm{z}}^{\max }$ is constrained to be positive using a meta-parameter ${\mathcal{F}}_{K_{\mathrm{z}}^{\max }}$ and the softplus function as $K_{\mathrm{z}}^{\max }=\log (\mathrm{1}+\exp ({\mathcal{F}}_{K_{\mathrm{z}}^{\max }}\left)\right)$. We ensure that $K_{\mathrm{z}}^{\min }\le K_{\mathrm{z}}^{\max }$ using a meta-parameter ${\mathcal{F}}_{K_{\mathrm{z}}^{\min }}$ as $K_{\mathrm{z}}^{\min }=K_{\mathrm{z}}^{\max }\mathit{\sigma }\left({\mathcal{F}}_{K_{\mathrm{z}}^{\min }}\right)$. We ensure that the shape parameter a is restricted between 0.5 and 1.5 using a meta-parameter ℱ_a as $a=\mathrm{0.5}+\mathrm{1.5}\mathit{\sigma }\left({\mathcal{F}}_a\right)$. The constrained super-elliptical function allows the calibration to identify a general nonlinear decrease in vertical hydraulic conductivity with increasing thaw depth.

Figure 2(a) Map; (b) hillshade relief with well locations; and (c) photo of the Barrow Environmental Observatory (BEO) Site A near Utqiaġvik, Alaska. Photo courtesy of Bob Busey.

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 low-centered 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) low-magnitude evapotranspiration signal, its effect on our calibration is relatively insignificant compared to the sporadic precipitation events that drive large-scale 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 pan-Arctic survey of ice-wedge polygon diameters does not to our knowledge currently exist, researchers have provided general characteristics based on extensive observations. Leffingwell (1915) states that ice-wedge polygons have an estimated average diameter of around 15 m. Liljedahl et al. (2016) indicate that low-centered polygons (including rims) can have diameters from 5–30 m. Abolt et al. (2018) state that polygonal formations in ice-wedge 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 low-centered polygon transition to high-centered polygons, at which time polygon-center 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 active-layer 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 active-layer thicknesses from 17–47 cm. Based on extensive ground-penetrating radar surveys near Utqiaġvik, Alaska, in 2013, Jafarov et al. (2017) document that average active-layer thickness was 41 cm with a standard deviation of 9 cm. Even for polygons with substantial active-layer thickness, it is important to understand the change in drainage throughout the season by evaluating early-season 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 ice-wedge polygons. Based on tracer arrival times, Wales et al. (2020) estimated horizontal conductivities of approximately 0.7–84 m d⁻¹ for a low-centered polygon and 0.01–0.3 m d⁻¹ for a high-centered polygon. Wales et al. (2020) stated that these are considered lower-bound estimates of horizontal conductivity since, based solely on breakthrough times, the approach is unable to isolate horizontal- and vertical-flow 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⁻¹, and vertical conductivities were around 0.2–0.4 m d⁻¹, 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 lower-bound 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⁻¹, 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 ice-wedge 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 higher-conductivity layers, and the effective vertical hydraulic conductivity captures flow in series across multiple layers and is dominated by the lower-conductivity 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 ice-wedge 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 non-existent 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⁻¹.

Figure 3(a) Progression of center water level calibration. The (1) base calibration, (2) calibration with a precipitation multiplier, and (3) calibration with vertical hydraulic conductivity as a function of thaw depth (K_z=f(D)). The measured polygon-center and trough water levels are plotted for reference. Precipitation is plotted on the right y axis. (b) Measured thaw depth and calibrated vertical hydraulic conductivity as a function of thaw depth.

Download

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 non-dimensional 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 non-dimensional 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 polygon-center 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 polygon-center 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 non-dimensional heads), and their effect on drainage timing (based on non-dimensional depletion curves) is straightforward, shifting the non-dimensional exponential drainage curve upwards or downwards. In other words, using non-dimensional variables is a powerful approach to gain intuition into the fundamentals of inundated ice-wedge polygon drainage irrespective of variable magnitude.

3.1 Calibration to field measurements

We present the best-fit polygon-center water levels for the calibration cases along with the measured values in Fig. 3. The two dominant drivers of the polygon-center 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²=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⁻¹; refer to Table 1) is extremely large and non-physical.

Considering that precipitation will run off from rims and collect in the polygon-center 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²=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⁻¹ 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²=0.96), as shown in Fig. 3 (red line). K_r is well within the estimated limits of 0.7 to 84 m d⁻¹ based on a tracer test at a nearby low-centered polygon conducted by Wales et al. (2020). Based on a rim width of 1 m, the discharge conductance of κ=3.3 d⁻¹ would correspond to a drainage boundary layer hydraulic conductivity of 3.3 m d⁻¹, which is also physically reasonable as opposed to the first two calibration estimates.

Table 1Calibrated parameters values, root mean square error (RMSE), coefficient of determination (R²) for calibration cases, and parameter standard errors for case (3). Note that the coefficient of determination is equivalent to the Nash–Sutcliffe efficiency here. Dashes indicate that the parameter was not part of the calibration.

Download Print Version | Download XLSX

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 meta-parameters) 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 polygon-center water level, and precipitation multiplier indicate that they are well constrained by the calibration.

The calibration verifies that the model is able to capture ice-wedge polygon drainage characteristics. In the next sections, we perform sensitivity analyses using non-dimensional forms of this verified analytical solution to gain insights into ice-wedge polygon drainage characteristics. The use of non-dimensional solution snapshots eliminates the need to consider the precipitation multiplier and thaw-depth-dependent vertical hydraulic conductivity explicitly. Instead, their effects are implicit in the relative differences between snapshots.

Figure 4Effect of polygon aspect ratio on polygon drainage with isotropic hydraulic conductivity. Plots along the left contain polygon radial transect head contours (colored lines) and filled stream tubes (gray regions) for several polygon aspect ratios (radius/thickness). The gray shaded region denotes the portion of the transect accessed by 95 % of the flow. The plots along the right contain corresponding gray rings indicating the surface area where 95 % of the polygon flow infiltrates. Each plot along the left contains a rectangle drawn to the actual proportions for the given polygon aspect ratio. In all cases, the anisotropy coefficient (horizontal conductivity/vertical conductivity) is fixed at unity.

Download

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 high-aspect-ratio 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.

Figure 5Effect of hydraulic conductivity anisotropy coefficient on polygon drainage. Plots along the left contain polygon radial transect head contours (colored lines) and filled stream tubes (gray regions) for several anisotropies (horizontal conductivity/vertical conductivity). The gray shaded region denotes the portion of the transect accessed by 95 % of the flow. The plots along the right contain corresponding gray rings indicating the surface area where 95 % of the polygon flow infiltrates. Each plot along the left contains a rectangle drawn to the actual proportions for the given polygon aspect ratio. In all cases, the polygon aspect ratio (radius/thickness) is fixed at 10.

Download

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 ice-wedge polygon drainage pathways.

Figure 6Effect of polygon aspect ratio on polygon drainage with high-hydraulic-conductivity anisotropy coefficient ($K_{\mathrm{r}}/K_{\mathrm{z}}=\mathrm{100}$). Plots along the left contain polygon radial transect head contours (colored lines) and filled stream tubes (gray regions) for several polygon aspect ratios (radius/thickness). The gray shaded region denotes the portion of the transect accessed by 95 % of the flow. The plots along the right contain corresponding gray rings indicating the surface area where 95 % of the polygon flow infiltrates. Each plot along the left contains a rectangle drawn to the actual proportions for the given polygon aspect ratio. In all cases, the anisotropy coefficient (horizontal conductivity/vertical conductivity) is fixed at 100.

Download

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.

Figure 7Contour map of the percent of the polygon access by 95 % of the flow as a function of polygon aspect ratio and anisotropy coefficient. Locations associated with other figures are indicated by labeled points. By choosing a constant value for R ($L=R/x$), horizontal transects parallel to the x axis represent the effect of thaw depth increasing from right to left, while they represent increasing polygon width from left to right for constant values of L (R=Lx). Note that κ=5 d⁻¹ in these calculations and therefore the location associated with Fig. 11b cannot be indicated.

Download

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).

Figure 8Depletion of ponded water height over time (H_c(t)) for (a, b) alternative polygon aspect ratios (radius/depth ($R/L$)) and (c) anisotropy coefficient ($K_{\mathrm{r}}/K_{\mathrm{z}}$). Plots (a) and (b) have an anisotropy coefficient fixed at unity, with plot (a) having high conductivity ($K_{\mathrm{r}}=K_{\mathrm{z}}=\mathrm{0.5}$ m s⁻¹) and plot (b) low conductivity ($K_{\mathrm{r}}=K_{\mathrm{z}}=\mathrm{0.05}$ m s⁻¹). The aspect ratio is fixed at 10 in plot (c). The curves also describe the depletion of ponded water volume over time (V(t)). Dashed lines indicate the characteristic times (t_L) in each case. Note that the orange line in (b) and (c) are identical, $R/L=\mathrm{10}$ and $K_{\mathrm{r}}/K_{\mathrm{z}}=\mathrm{1}$.

Download

3.3 Ponded water depletion

Depletion curves for the non-dimensional 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 y-axis 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.

Figure 9Characteristic time as a function of (a) polygon aspect ratio and (b) anisotropy coefficient. The $K_{\mathrm{r}}/K_{\mathrm{z}}=\mathrm{1}$ line in plot (a) corresponds to characteristic times in Fig. 8b. Note that the trends in (a) are nearly linear and on naturally scaled axes, while the trends in (b) are also nearly linear on log-transformed axes.

Download

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.

Figure 10Contour map of the characteristic time of polygon drainage as a function of polygon aspect ratio and anisotropy coefficient. Locations associated with other figures are indicated by labeled points. By choosing a constant value for R ($L=R/x$), horizontal transects parallel to the x axis represent the effect of thaw depth increasing from right to left, while they represent increasing polygon width from left to right for constant values of L (R=Lx). Note that κ=5 d⁻¹ in these calculations, and therefore, the location associated with Fig. 11b cannot be indicated.

Download

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.

Figure 11Attaining the same mathematical solution by modifying polygon aspect ratio ($R/L$), anisotropy coefficient ($K_{\mathrm{r}}/K_{\mathrm{z}}$), and rim conductance (κ). The left plots of (a) and (b) contain polygon radial transect head contours (colored lines) and filled stream tubes (gray regions) for several polygon aspect ratios (radius/thickness). The gray shaded region denotes the portion of the transect accessed by 95 % of the flow. The right plots of (a) and (b) contain corresponding gray rings indicating the surface area where 95 % of the polygon flow infiltrates. Plots (a) and (b) contain a rectangle drawn to the actual proportions for the given polygon aspect ratio. The ponded water height and volume depletion curves for both cases are presented in (c).

Download

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 non-dimensional 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 non-modified 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 outlet-boundary-interface-limited 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_r1K_z1=K_r2K_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², 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², 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.

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 real-world applicability of the sensitivity analysis based on comparisons of steady-state 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 thaw-depth-dependent 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 sub-centimeter RMSE of 0.37 cm and an R² (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 thaw-depth-dependent vertical hydraulic conductivity in the model. The other parameters (discharge conductance, initial polygon-center 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 ice-wedge 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 intra-polygon 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 active-layer depth are relatively common in polygonal terrain. It is unlikely that these sometimes large geochemical-depth gradients would be preserved with significant amounts of deep flushing. Inter-annual 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 ice-wedge 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 ice-wedge top thawing, which would promote low- to high-centered 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 ice-wedge top as well (Fig. 6d). Similarly, increasing seasonal thaw depth and inter-annual active-layer thickness will lead to polygons with less focused advective heat transport towards ice-wedge tops, potentially providing a negative feedback on ice-wedge degradation. These results provide an additional perspective on ice-wedge 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 inter-annual active-layer 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, large-scale 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 ice-wedge 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 ice-wedge 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 ice-wedge polygon geometry will best approximate drainage in nearly symmetrical hexagonal ice-wedge polygons, while non-symmetric or square ice-wedge 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 non-idealized 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 ice-wedge polygon drainage and timing.

As the analytical solution applies to ponded conditions, it does not apply to freeze-up 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 ice-wedge 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 freeze-up (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 non-inundated low-centered polygons. Since it is based on having a ponded center, the model is also not applicable to high-centered 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 ice-wedge 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, Cresto-Aleina et al. (2013) analyze low-centered 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 ice-wedge 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 process-rich 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 process-rich, complex 1D model to calibrate and perform uncertainty and sensitivity analyses of thermal-hydrology models of ice-wedge 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 ice-wedge thaw on individual polygons using a full-physics thermal-hydrology 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 full-physics thermal-hydrology model of ice-wedge 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 ice-wedge polygon drainage, and provides new modeling and experimental research directions based on the idealized base cases considered here.