New insights into the drainage of inundated Arctic polygonal tundra using fundamental hydrologic principles

The pathways and timing of drainage from inundated ice-wedge polygon centers in a warming climate have important implications for carbon flushing, advective heat transport, and transitions from carbon dioxide to methane dominated emissions. This research helps to understand this process by providing the first in-depth analysis of drainage from a single polygon based on fundamental hydrogeological principles. We use a recently developed analytical solution to evaluate the effects of polygon aspect ratios (radius to thawed depth) and hydraulic conductivity anisotropy (horizontal to vertical hydraulic 5 conductivity) on drainage pathways and temporal depletion of ponded water heights of inundated ice-wedge polygon centers. By varying the polygon aspect ratio, we evaluate the effect of polygon size (width), inter-annual increases in active layer thickness, and seasonal increases in thaw depth on drainage. One of the primary insights from the model is that most inundated ice-wedge polygon drainage occurs along an annular region of the polygon center near the rims. This implies that inundated polygons are most intensely flushed by drainage in an annular region along their horizontal periphery, with implications for 10 transport of nutrients (such as dissolved organic carbon) and advection of heat towards ice wedge tops. The model indicates that polygons with large aspect ratios and high anisotropy will have the most distributed drainage. Polygons with large aspect ratio and low anisotropy will have their drainage most focused near the their periphery and will drain most slowly. Polygons with small aspect ratio and high anisotropy will drain most quickly.


Introduction
Polygonal tundra occurs in continuous permafrost landscapes lacking exposed bedrock or active sedimentation (Brown et al., 1997;MacKay, 2000) and is estimated to occupy around 250,000 km 2 , or around 3% of the Arctic landmass predominantly along the north coasts of Alaska and Siberia (Minke et al., 2007). It is rich in organic carbon (Tarnocai et al., 2009;Hugelius et al., 2014) which is poised to release to the atmosphere as carbon dioxide or methane under a warming climate causing a 20 potential feedback to climate change (Schuur et al., 2008). The microtopography of polygonal tundra collects precipitation and snow melt events resulting in a slow release from the landscape through the thawed subsurface to surface water drainage Whole polygon represents an idealized pie wedge slice of a low-centered polygon including a wedge of the ponded center, rim, and trough. Equipotential hydraulic head lines are denoted as h(r, z) and stream lines as Ψ(r, z). L is the depth of the thawed soil layer and R is the radius of the polygon center. The ponded water height (Hc(t)) and trough water level (Ht) are noted.
networks. This process has important implications for the release of dissolved organic matter from polygon tundra landscapes flushed from the subsurface to surface waters and transport of heat through the subsurface. The role that polygonal tundra regions play in global terrestrial biogeochemical feedbacks as the Arctic warms necessitates understanding their hydrologic 25 flow regimes, including carbon fluxes from the landscape.
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 re-cracking. Over many cycles, this process leads to the growth of subsurface ice wedges connected in polygonal patterns known as ice-wedge polygons . Ice-wedge polygons vary in size from several meters to a few tens of meters in diameter and develop 30 over time frames of hundreds to thousands of years (Leffingwell, 1915;Lachenbruch, 1962;MacKay, 2000;Abolt et al., 2018).
Thermal expansion of the ice and soil during summer months often results in warping the soil strata forming parallel rims on both sides of the ice wedge and a trough directly above the ice wedge (refer to Figure 1). The peat and mineral soil layering often approximately follow the surface topography, resulting in raised mineral soil under the rims which can impede drainage due to the mineral soils relatively lower hydraulic conductivity (Hinzman et al., 1991). Ice-wedge polygons with well-formed rims 35 and troughs with a distinct, central topographic depression are referred to as low-centered polygons. Ponding occurs frequently in the centers because the rims can trap water during snow melt 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). As observed by Helbig et al. Peterson, 1980;Jorgenson et al., 2015;Wolter et al., 2016), and (4) ground surface deformation (Mackay, 1990;MacKay, 2000;Raynolds et al., 2014;Oehme, 2019).
The microtopographic features of polygonal tundra result in sharply contrasting spatial thermal and hydrologic conditions (Boike et al., 2008;Zona et al., 2011;Helbig et al., 2013) leading to sharply contrasting spatial biogeochemical conditions Norby et al., 2019). Standing water in polygonal centers show distinctly different chemistry compared 60 to surface water in troughs and ponds (Zona et al., 2011;Koch et al., 2014Koch et al., , 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 inundated center drainage to troughs. The importance of representing carbon respiration response to changing hydrologic conditions in arctic landscape has motivated the development of complex biogeochemical modules in Earth system models that are driven by hydrologic conditions (Grant et al., 2017;Bisht The geometric shape of low-centered polygons along with soil hydraulic properties (for example, hydraulic conductivity) affect the distribution of hydraulic heads that control the pathways and timing of inundated ice-wedge polygon drainage. In 75 this paper, we use fundamental hydrogeological principles (Harr, 1962;Cedergren, 1968;Freeze and Cherry, 1979;Bear, 1979) to understand effects of geometry and hydraulic properties on inundated ice-wedge polygon drainage. We investigate the pathways and timing of inundated polygon drainage using a 3D-axisymmetric analytical solution of groundwater heads. This approach idealizes the thawed subsurface of an inundated low-centered 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; annular 80 ring defined by a line from (R, 0) to (R, L) in Figure 1). The depth (L) and radius (R) of the thawed soil layer of the polygon center and horizontal (radial) and vertical hydraulic conductivity are variable within the solution. It is important to note that we are focusing on drainage of inundated low-centered polygons here, and that biogeochemical investigations suggest that the drainage of non-inundated low-centered, transition, and high-centered polygons may have different characteristics Newman et al., 2015;Wainwright et al., 2015;Wales et al., 2020).

85
This model provides fundamental insights into the pathways and timing of inundated ice-wedge polygon drainage for polygons with different geometries and degrees of anisotropy. Varying the geometry of the cylinder not only allows us to capture differences in drainage between different sized polygons (in other words, polygons with different radii), but also seasonal and inter-annual variations in polygon thaw-layer depth (in other words, as the thawed soil layer increases during the summer or as the active layer increases from year to year). By allowing for different hydraulic conductivities between the horizontal and 90 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 high conductivity peat layers overlying lower conductivity mineral soil results in horizontal watershed drainage (McDonnell et al., 1991;Brown et al., 1999;Quinton and Marsh, 1999). Helbig et al. (2013) found that in polygon tundra, the peat layer quickly redistributed precipitation events 95 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. While the cylindrical geometry and anisotropy considered here do not cover the potential variations that are known to be present in ice-wedge polygons (for example, anomalous heterogeneities), the impact of those variations will cause deviations around the base case scenarios considered here.
Although there have been numerous field observations of inundated low-centered polygon drainage to troughs (Helbig et al.,100 2013; Koch et al., 2014;Wales et al., 2020), the research here is, to our knowledge, the first in-depth investigation into the pathways and timing associated with inundated ice-wedge polygon drainage based on fundamental hydrologic principles. Aside from the field tracer experiments of Wales et al. (2020), previous research has primarily focused on 1D vertical hydrothermal effects within ice-wedge polygons (Atchley et al., 2015;Harp et al., 2015;Atchley et al., 2016) or ice-wedge surface hydrology (Boike et al., 2008;Jan et al., 2018a, b). Other researchers have investigated the hydrology of multiple polygons without 105 investigating the drainage pathways within ice-wedge polygons themselves (Nitzbon et al., 2019). Our analysis provides a new perspective on inundated polygon hydrogeology indicating the implications of polygon geometry and hydraulic conductivity anisotropy on drainage pathways and timing.

Model parameterization
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 currently exist to our knowledge, 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.  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 115 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 near 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 120 result in low-centered polygon transition to high-centered polygons, at which time center inundation will cease to occur due to rim collapse. Therefore, there is a limit to the thaw depth of interest. 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 to 2009 had active layer thicknesses from 17 to 47 cm. Based on extensive ground penetrating radar surveys near Utqiaġvik, Alaska in 125 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 can be captured by larger aspect ratios.
Based on these considerations, we evaluate aspect ratios from 2.5 to 20. For example, given a thawed 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 130 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 thawed depth of 0.5 m, an aspect ratio of 20 would represent a polygon with a 10 m radius center.
While our selected aspect ratios and anisotropies may not cover all possible scenarios, they do cover a broad enough range of scenarios to draw insights into their effect on inundated ice-wedge polygon drainage.
Comprehensive measurements of hydraulic conductivity anisotropy are lacking from ice-wedge polygons. Based on tracer 135 arrival times, Wales et al. (2020) estimated horizontal conductivities of approximately 0.7 to 84 m/day for a low-centered polygon and 0.01 to 0.3 m/day 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) perform laboratory measurements of anisotropy on small samples from peat soils where horizontal conductivities were around 6-7 m/day and vertical conductivities were around 0.2-0.4 m/day 140 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 freeze-thaw dynamics that current ice-wedge polygons will have been subjected to in recent times, 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/day, primarily focusing on anisotropic values from 0.1 145 to 100. These values allow us to investigate the effects of anisotropy on the drainage pathways and timing using hydraulic conductivities consistent with currently available measurements.

Model overview
We use recently developed 3D-axisymmetric analytical solutions derived in Zlotnik et al. (2020), and briefly described in Appendix A, of hydraulic heads and the stream function in the thawed soil layer below a polygon center (dashed rectangle 150 in Figure 1) to investigate inundated low-centered polygon drainage pathways. In order to generalize the results with regard to polygon center and trough ponded heights, we use non-dimensional forms of the hydraulic head and stream function (refer to equation A6 and A7). 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 absolute magnitude of these patterns as the drainage process proceeds. In order to evaluate 155 drainage pathways, we plot flow nets, as illustrated in Figure 1, composed of lines of equal non-dimensional hydraulic head (h * (r * , z * )), referred to as equipotentials, and contours of the stream function (Ψ * (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 by 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 160 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 Figure 1).
We normalize the stream function Ψ(r * , z * ) from 0 to 1, and similar to heads, plot stream function contours (streamlines) 165 from 0.05 to 0.95 by increments of 0.1. Given this set of streamlines, 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), convey 5% of the drainage each.
In order to quantitatively evaluate the spread of drainage within the polygon, 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 170 time can be expressed by a simple exponential decay function as where t L is the characteristic time of drainage (refer to equation A10), defined as the time when the ponded height is 1/e times, or ∼37% of, its original height H c,0 . The dimensional ponded height H c (t) can be obtained as 175

Results
Here, we present drainage flow nets for various aspect ratios and anisotropies, polygon center 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. We compare the effect of aspect ratio on drainage pathways in an isotropic and a highly anisotropic (K r /K z = 100) polygon. The drainage pathways are also evaluated for various anisotropies holding the aspect

Drainage pathways
The effect of aspect ratio on drainage pathways when the hydraulic conductivities are isotropic (K r /K z = 1) is shown in Figure 2. 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 190 words, 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, 95% of the drainage is focused within around 5% of the polygon volume, while for an aspect ratio of 2.5, the drainage is spread over around 29% 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 Figure 2 indicate that 195 throughout the thaw season, as the thaw depth increases, the drainage path will spread out towards the middle of the polygon center. Similarly, Figure 2 can be used to evaluate drainage pathways of polygons of different widths but similar thaw depths.
In this context, Figure 2 indicates that wider polygons will have more focused drainage, while drainage for the small polygons will be more dispersed.
The effect of anisotropy on drainage pathways when the aspect ratio is held constant (R/L = 10) is shown in Figure 3.

200
As anisotropy increases (in other words, the more that horizontal conductivity dominates), the region accessed by drainage flow becomes larger. Similar to a decreasing aspect ratio in Figure 2, increasing anisotropy leads to a larger radial extent of the accessed region, while the vertical extent is nearly unaffected. When the vertical conductivity is ten times the horizontal (top plot), only around 3% of the polygon volume is accessed by 95% of the drainage. If horizontal conductivity is 100 times   The effect of aspect ratio on drainage pathways when the hydraulic conductivities are highly anisotropic (K r /K z = 100) is shown in Figure 4. In this case, contrary to the isotropic case in Figure 2, as the aspect ratio decreases, the accessed volume generally decreases. However, the largest accessed volume is for an aspect ratio of 10, not 20. This nuance in the dependence of    of the accessed volume as aspect ratio decreases. By comparing Figures 2 and 4, it is also apparent that the volume accessed by 95% of the drainage is generally larger with higher anisotropy.
To gain a global perspective on the trends in drainage pathways with aspect ratio and anisotropy, Figure 5 maps the percent of the polygon accessed by 95% of the drainage as a function of aspect ratio and anisotropy. This approach illustrates the overall structure of the combined effect of aspect ratio and anisotropy on accessed volume (focused versus dispersed drainage).

215
The trend in accessed volume with respect to aspect ratio in Figure 2 is represented along the anisotropy=1 transect in Figure 5, while the trend in Figure 4 is represented by the anisotropy=100 transect. There is a spreading ridge-like structure in the accessed volume with increasing aspect ratio and anisotropy with decreasing accessed volume on either side of this curved ridge line. The ridge represents the optimal balance of radial and vertical extension of the accessed volume region, as described in the previous paragraph.

220
This structure explains the increasing accessed volume with decreasing aspect ratio in Figure 2 and the maximum accessed volume at aspect ratio of 10 in Figure 4. Similarly, the trend in accessed volume with respect to anisotropy is represented by Ponded volume the aspect ratio=10 transect in Figure 3. The drainage flow is most spread out (largest accessed volume) when the aspect ratio is large and the anisotropy is high (yellow region in the upper left of Figure 5). The drainage is the most focused (least accessed volume) when the aspect ratio is large and the anisotropy is low (dark blue region in the lower left of Figure 5).

Ponded water depletion
Depletion curves for the non-dimensional ponded water height in the polygon center for various aspect ratios and anisotropies are shown in Figure 6. Note that since the depletion of the volume of ponded water is directly related to the ponded water height, the plots in Figure 6 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 1/e ≈ 0.37 its initial 230 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 higher anisotropy will drain faster than those with low anisotropy.
It is also apparent that the increase in drainage time with aspect ratio is nearly linear, while for anisotropy, it is exponential. This is further illustrated in Figure 7, where we plot the trend in characteristic times as a function of aspect ratio for various 235 anisotropies (top plot) and as a function of anisotropies for various aspect ratios (bottom plot). It is apparent that 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 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 240 anisotropy is shown in Figure 8. The fastest (shortest) drainage times are achieved with low aspect ratio and high anisotropy, while the slowest (longest) drainage times are achieved with high aspect ratio and low anisotropy. 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 ratio, those with preferential horizontal flow will drain most quickly, while those with preferential vertical flow will drain most slowly.

Counteracting effect 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, anisotropy stretches the domain by multiplying the non-dimensional radius by (K z /K r ) (refer to definition of r * in equation 1). As a result, the effect of increasing anisotropy 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 250 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 identical effect as dividing the anisotropy by four. However, this is not the case because the solution involves a Biot parameter (Bi) which defines a 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 resistance to drainage 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, one would also need to ensure that Bi is not modified in the process.

260
However, given two solutions, this requires satisfying the conflicting conditions that K z1 /K r1 = K z2 /K r2 and K r1 K z1 = K r2 K z2 (refer to equations 1 and 4), 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 to achieve an identical mathematical solution of drainage. An example where identical mathematical solutions for drainage are obtained is provided in Figure 9, where the bottom plot is obtained by taking the properties of the top plot and dividing the aspect ratio by 8, dividing the 265 anisotropy 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 4L in the top case and around 0.5L in the bottom case).
In the bottom plot in Figure 9, we demonstrate that despite these differences, the depletion curves are identical for both cases.

270
This result is due to our choice to modify R to change the aspect ratio and modify K r to change the anisotropy. If either L or K z are chosen to modify the aspect ratio or anisotropy, respectively, the depletion curves would not necessarily by 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 as above), the drainage flow net would still be mathematical equivalent, but the characteristic time would be 8 times shorter than in the original case (refer to equation A10). Therefore, while it is possible to obtain mathematically equivalent 275 drainage patterns with counteracting modifications to aspect ratio, anisotropy, and outflow resistance, the temporal depletion will only be equivalent if R and K r are used to modify the aspect ratio and anisotropy, respectively.

Discussion
Our analysis provides new insights into the manner in which geometry and anisotropy effect how inundated ice-wedge polygons retain and slowly release water from their centers to their troughs, which form the drainage network of polygonal tundra 280 landscapes. Using a mathematical representation of inundated ice-wedge 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), we quantify the sensitivity of inundated ice-wedge polygon drainage to representative polygon sizes, inter-and intra-annual changes in thaw depth, and preferential flow (hydraulic conductivity anisotropy).
A key result of this study is that the geometry and anisotropy of the polygon subsurface have a significant effect on the region 285 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 aspect ratio and anisotropy. In addition, the majority of the ponded water will flow to an outer ring of the polygon before infiltrating into the subsurface. Field observations have 290 indicated not only the existence of intra-polygon biogeochemical diversity (Zona et al., 2011;Newman et al., 2015), but 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 Newman et al., 2015;Wales et al., 2020;Plaza et al., 2019) effecting dissolved organic matter and mineral flushing from the subsurface of polygon centers to the surface waters of troughs. This is important not only simply 295 with regard to 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 300 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; Figure 2 and 4). However, with high anisotropy, 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 contract towards the ground surface ( Figure 4). 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 low aspect ratio and high anisotropy, the dissolved 305 organic matter and nutrients closer to the surface will be flushed more than the deeper stores of dissolved organic matter and nutrients. 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 high anisotropy. This research reinforces the need for field studies on anisotropy and preferential flow in polygon landscapes in order to better understand the hydrologic transitions and feedbacks that will occur in a warming climate.

310
The drainage pathways and timing presented here are based on fundamental hydrogeological principles (Harr, 1962;Cedergren, 1968;Freeze and Cherry, 1979;Bear, 1979) using a generalization of ice-wedge polygon geometry and hydraulic properties that allow fundamental insights of generalized drainage pathways and timing of inundated ice-wedge polygons to be obtained. Our solutions capture the fundamental physical forces that ponded water exerts on a cylindrical porous disc with drainage allowed radially through its sides. While the cylindrical geometry and anisotropy considered here do not cover all 315 potential variations present in ice-wedge polygons, the impact of those variations will cause deviations around the base case, idealized scenarios considered here. For example, anomalous heterogeneities will warp the base case hydraulic head equipotentials and stream lines. The cylindrical idealization of ice-wedge polygon geometry will best approximate drainage in nearly symmetrical hexagonal ice-wedge polygons, while non-symmetric, square ice-wedge polygon drainage will deviate most from that shown here. Ice-wedge polygons with focused, deep annular flow towards their periphery may also direct warmer surface waters towards the top of ice-wedges, resulting in enhanced ice-wedge degradation (Wright et al., 2009). Based on our calculated drainage pathways, advective heat transport will be most pronounced in large aspect ratio polygons with low anisotropy (refer to the top plots in Figures 2 and 3). 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 high anisotropy as in the upper plots in Figure 4). However, small polygons with high anisotropy will restrict the drainage flow near the ground surface, potentially reducing the advective transport of heat to the ice-wedge top as well (bottom plot in Figure 4).
Small polygons with deeply thawed soil layers (low aspect ratios) and high horizontal preferential flow (high anisotropy) 335 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 would 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. Of course, for a given location, other factors such as regional flow patterns, large scale topography, etc. will influence the regions overall drainage timing. However, along with these other factors, our 340 analysis indicates that aspect ratio will have a nearly linear positive relationship while anisotropy will have an exponential negative relationship with drainage timing.

Conclusions
-The majority of drainage from inundated ice-wedge polygon centers occurs along an annular region along their radial periphery; however, polygon aspect ratio and hydraulic conductivity anisotropy significantly impact the drainage path-345 ways.
-A combination of high aspect ratios (wide, shallow polygons) and high anisotropy (preferential horizontal flow) result in the greatest spreading of drainage flow and largest fraction of the polygon volume being accessed by drainage flow.
-A combination of high aspect ratios (wide, shallow polygons) and low anisotropy (preferential vertical flow) result 350 in the greatest focusing of drainage flow and smallest fraction of the polygon volume being accessed by drainage flow.
-Combinations of aspect ratio and anisotropy have counteracting effects of radial versus vertical extension/contraction of drainage pathways, producing non-monotonic relationships between aspect ratio/anisotropy and accessed volume (ridgeline in accessed volume response surface in Figure 5).

355
-The characteristic time for a polygon to drain has an approximately positive linear relationship with aspect ratio; in other words, wide, shallow polygons drain slowly while small, deep polygons drain quickly.
-The characteristic time for a polygon to drain has a negative exponential relationship with anisotropy; in other words, preferential horizontal flow leads to exponentially faster drainage.
Code availability. Matlab scripts of the analytical solutions are provided as supplementary information.
We idealize the subsurface below the center region of a low-centered polygon as a thin cylinder with radius in the horizontal direction and length in the vertical direction (refer to Figure 1). Based on this approximation, the initial (steady-state) hydraulic heads (h(r, z, 0)) in an inundated polygon will satisfy the following equation: where h is the hydraulic head, r is the radial coordinate, z is the depth coordinate (positive in the downward vertical direction),

370
K r is the horizontal (radial) hydraulic conductivity, and K z is the vertical hydraulic conductivity. The system is fully specified by ensuring that (1)

380
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 resistance to flow across the drainage interface to the trough (in our case, the hydraulic resistance of the soil layer under the rim).
The solutions can be verified by direct substitution of equations A6 and A7 into the boundary value problem defined by equations A1 to A5.

395
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 c (t)/H c,0 = 1/e; in other words, when the ponded height is approximately 37% of its initial. The 400 function F (Bi, R * ) can be evaluated as F (Bi, R * ) = R * 0 ∂h * (r * , 0) ∂z * r * dr * = 2 ∞ n=1 tanh(λ n ) λ n (λ n /Bi) 2 + 1 .
The depletion curve can be expressed in non-dimensional terms as where non-dimensional time t * = t/t L . The solution can be verified by direct substitution of equation A9 into equation A4.