Thermal contraction crack polygons in Nunavik (northern Quebec): Distribution and development of polygonal patterned ground

We evaluated the spatial distribution and morphological variability of thermal contraction crack polygon (TCCP) networks across Nunavik, a 440,000‐km2 region of northern Québec that spans the northward transition from discontinuous to continuous permafrost. A population of 4,567 TCCP sites was sampled and analyzed from 80,737 georeferenced high‐resolution aerial photographs and 264,504 km2 of ESRI satellite basemaps. For each site, six parameters were inventoried and compiled into a database: (a) network geometric arrangement; (b) intersection angles; (c) number of subdivisions and nested polygons (referred to as generations of development); (d) dominant polygon morphology; (e) surficial geology; and (f) vegetation cover. Statistical analyses of the tabulated data revealed a strong association between Holocene glacial, glacio‐fluvial, fluvial, marine, and organic landforms and the different intersections angles in the networks, providing insight into how the processes of thermal contraction cracking function and manifest geomorphically across varied permafrost landscapes. Orthogonal polygons (intersection angle of 90°) dominate on flat terrains where the thermo‐mechanical stresses are probably spatially homogeneous. Hexagonal (angles of 120°) and poorly structured polygons tend to form where topography variability probably generates heterogeneous heat flow patterns and thermo‐mechanical stresses in the ground, resulting in irregular cracking patterns.


| INTRODUCTION
Tundra polygons are the most widespread periglacial landforms in subarctic and arctic environments, and the ice wedges that underlie them account for an important fraction of ground ice in permafrost. 1,2 Polygons in the periglacial environment develop due to thermal contraction cracking under low winter temperatures and during spells of rapid cooling that generate vertical frost cracks that may reach depths of up to 10 m. 3 Frost cracks organize themselves in various geometric patterns (e.g., quadrangles [90 , T-junction] or hexagons [120 , Yjunction]) 4 depending on the geometry of the stress field in the terrain, following theoretical thermo-mechanical principles. 5 Recent observations and field measurements with instrumentation using breaking electrical wires, extensometers, accelerometers, thermistors, and data loggers have confirmed the theory and yielded some empirical information on the temperature regimes that regulate frost cracking. [6][7][8][9][10][11][12][13][14][15] Contraction cracks typically infill with snowmelt water in spring, and refreezing of this snowmelt leads to the formation of ice wedges. 6 However, in some environments, depending on local conditions such as soil moisture and sediment availability, the cracks may be filled with sand, organic matter, mud, or mixtures of soil and ice, feeding the formation of sand wedges, composite wedges (mixed), soil wedges, and soil veins. 1 Because the methodology in this study is based mostly on surface pattern classification and mapping over a vast territory with a limited number of field observations, we use herein the generic term "thermal contraction crack polygon" (TCCP) as it covers the different varieties of documented crack-filling polygons.
Nunavik (northern Quebec) is a region with a strong climate gradient transitioning between boreal forest and tundra. Moreover, the region has a complex Holocene history with soil materials derived from various surficial geological settings and with terrain exposure of various ages associated with deglaciation, post-glacial land emergence, Holocene climate changes, and vegetation history (e.g., tree-line fluctuations, forest fires and recent shrubification).
Climate warming and subsequent permafrost thaw pose a significant risk to infrastructure in Arctic communities. 16,17 Environments with TCCP networks are considered particularly sensitive to climatic disturbances and thawing permafrost and melting ice wedges can create ponds, settlement, and thermo-erosional gullies, resulting in possible damage to roads, houses, and airport runways. [18][19][20][21] Despite abundant research in Nunavik on permafrost, TCCPs and on ground thermal contraction processes, no systematic mapping of TCCPs and analysis of their geomorphic and climatic controls has been undertaken before. As a result, the geological and ecological terrain conditions under which TCCPs formed and have evolved, as well as how they manifest across the landscape, are largely unknown. As most of the permafrost in Nunavik is epigenetic and aggraded originally in pre-existing Quaternary landforms, we can infer that the pre-existing geomorphological contexts exerted some physical control on the formation and the resulting morphologies of the TCCPs.
We propose that by classifying morphological and ecological differences at multiple polygon sites across a range of geological and climatic contexts, conclusions can be drawn about past and current cracking conditions in Nunavik, with some inferences on permafrost history and future terrain susceptibility to change. More specifically, the specific objectives aim to: 1. Classify the TCCPs in Nunavik based on their morphology.

Test statistical relationships between the morphologies of TCCPs
and Quaternary deposits and landforms. 3. Improve understanding of the spatial distribution of TCCPs and their morphological characteristics.
These data and analyses will provide a basis for deriving thermomechanical processes of TCCP formation under a range of common/ typical site conditions and inferring the impacts of Holocene ecological changes on TCCP activity and morphology.

| BACKGROUND
Frost contraction in permafrost terrain and opening of linear cracks forming polygonal patterns occur after freeze-back of the active layer, usually in midwinter when ground temperatures are below at least À10 C for prolonged periods and a ground cooling rate ranged between À0.1 and À0.6 C d À1 for two or more successive days ( Figure 1). 1,5,[10][11][12][13][14][15]22 Regardless of the type of crack infilling (i.e., ice, sand, organic matter), the networks share common configurations: a pattern of cracks forming often orthogonal (90 ) and, sometimes, hexagonal (120 ) surface patterns, which may be poorly organized or oriented with variable sizes of polygons. 5 The surface factors governing the formation of TCCPs (air and ground temperatures, vegetation cover, snow cover, and soil materials and moisture) have been abundantly discussed by previous researchers. 1,[5][6][7][8][9][10][11][12][13][14][15]22 For instance, Mackay 23,24 highlighted the relationship between air temperature and cracking in areas with low snow cover. In fact, heavy snowfall can inhibit cracking. Along the western Arctic coast, no ice new veins above the ice wedges were observed when the average snow depth exceeded about 60 cm.
Orthogonal patterns develop in an evolutionary sequence in which the primary cracks are intersected by secondary cracks that subdivide the stress field, thus forming right-angled, orthogonal patterns. 1,5,25 Primary cracks develop perpendicular to the general orientation of the deposit (e.g., a raised beach or a river terrace). Secondary cracks (often smaller) intersect them, cutting out quadrilaterals. These observations apply particularly well to newly emerged areas where it is possible to see the propagation and expansion of cracks in sediments ( Figure S1). For hexagonal networks (120 ), Lachenbruch 5 and Plug and Werner 26,27 suggest that the cracks develop simultaneously, probably from random points in a homogeneous soil with uniform thermo-mechanical properties.
Ice-wedge polygons are the most abundant, observed, and studied TCCPs. 1,6,[9][10][11][12][13]21 They are located exclusively in permafrost environments and cover almost a million square miles (2,600,000 km 2 ) in tundra and boreal forest across the Northern Hemisphere. 28 They occur widely in nonbedrock permafrost environments, such as lowland areas, tundra wetlands, alluvial plains, along hilltops, and on glacial deposits, 6,9,13,[29][30][31] and can have different morphologies, being either low-centered, flat, or high-centered. 9 Baydjarakhs, 32 also spelled baydzherakh, 20,33 are patterned fields of mounds resulting from advanced ice-wedge melting in high-centered polygons (see section 8 Terminological note). Sand wedges are known to develop in cold and arid polar regions, where there is little or no water supply and a very thin snow cover. 34 Sand wedges occur where mean annual air temperatures are between À4 and À8 C and mean annual precipitation is <100 mm. Black 35  to permafrost zones, but can also be found in seasonally frozen regions and are formed principally near the surface, in the active layer and the transition layer and do not propagate into permafrost. 1,[35][36][37][38][39][40][41] In Nunavik, soil wedges are in areas where mean annual ground temperatures are between 0 and À5 C. 38 3 | NUNAVIK Nunavik is the Inuit territory of northern Quebec, above the 55th parallel and extending northward to the Hudson Strait, with an area of 443,685 km 2 . The spatial distribution of permafrost in Nunavik has been classified into five zones: continuous (>90% of the ground surface), extensive discontinuous (50-90% of the ground surface), discontinuous and dispersed (10-50% of the ground surface), sporadic (1-10% of the land area), and isolated patches or residual permafrost (<1% of the land area) ( Figure 2). 38,42,43 In the Ungava Peninsula, permafrost is continuous and reaches depths of up to 630 m. 44 In southern Nunavik, permafrost occurs sporadically or in isolated patches, and is often limited to $10-20 m in thickness. Inland, bedrock geology consists of expanses of exposed Precambrian intrusive and metamorphic rocks, with surficial materials consisting of shallow (<2 m) and thick till covers and fields of glacial and glaciofluvial deposits. In the coastal regions that were covered by The post-glacial climate over the territory above the treeline in parts of Nunavik that was never submerged by post-glacial seas was conducive to the formation and the sustained presence of permafrost in what is now the continuous permafrost zone. However, the discontinuous permafrost parts of the region with shrub and forest tundra have a more complex ecological and climatic history that involves alternating periods of growth and decay of permafrost during the Late Holocene period and the 20th century. 10,38 The timing of these paleoclimate shifts were also corroborated through the stratigraphic interpretation of past ice-wedge activity near Salluit (northern Quebec, Nunavik), 10,45 and later backed by climate reconstruction by reverse modeling of deep temperature profiles at Raglan Mine. 44 In all cases, the Little Ice Age (LIA) (roughly dated from $1450 to 1900 AD) stands out as the probable coldest Late Holocene period in the region. 10 (Table 1).
When a polygonal network fulfilling the criteria established in Table 1

| Acquisition of aerial photographs (TCCP mapping)
We used GPS-georeferenced high-resolution oblique aerial photo-  coverage. The analysis of ESRI satellite imagery (1:10,000) enabled us to cover 264,504 km 2 of territory, representing more than half of the area of Nunavik ( Figure 3c). The resolution of the basemaps ranged from 0.10 to 1 m, which was sufficient to distinguish the parameters associated with the polygons. Each site with TCCPs was systematically analyzed according to the list of parameters in Table 1.
TCCPs located using the oblique aerial photographs were repositioned to their appropriate locations to avoid duplication with the TCCPs located via the ESRI basemaps. Area measurements for polygon networks were performed using the raster basemaps and the editor tool in ArcMap 10.7. The areas were calculated based on the Quebec Lambert Conformal Conic projection. To assess the influence of local and regional topography on TCCP parameters, topography and slope data were derived from the ArcticDEM, which is a 2-m digital elevation model. 52  To minimize the uncertainty of the mapping procedure, interpretation of the aerial photographs was done twice, once at the beginning of the project, and a year later before starting the statistical tests.

| Chi-squared test (χ 2 )
A statistical analysis of the parameters of the mapped polygons was applied to test whether associations existed between the different TCCP morphological parameters and between the TCCPs and the environmental factors. The association between two categorical variables can be assessed using the chi-square test (χ 2 ) (Equation 1). [53][54][55][56][57] The general formula for χ 2 is: where O is the observed cell frequency, E is the expected cell frequency, i runs from 1, 2, … to n, where n is the number of cells in the contingency table, and Σ stands for the sum of the results combined. 56 For instance, the χ 2 test was performed using the intersection angle and surficial deposit variables to determine whether these variables are associated (Table 2). This allows examination of whether the difference between the observed and expected frequencies (Equation 2) of two categorical variables (here, the variables are the network parameters) is due to chance. 53 In Table 2, the number of TCCPs forming in fluvial deposits with mixed intersection angles (observed: O) is 36; to find the expected frequency (E) for the χ 2 test, the following formula was applied: where Row marginal total is the sum of the row entries (445), Column marginal total is the sum of the column entries (359), and Overall total is the total sum of the sample (4,307). 56 Here, the expected value is 37.09. Once the values of all the expected frequencies are determined, we can deduce the χ 2 values. In this example case, the χ 2 value for The null hypothesis is rejected when the probability or P-value is less than the alpha level (P ≤ 0.05 or P ≤ 0.01), also called a significance level, chosen according to the degrees of freedom (df ) (Equation 3).
where R is the number of rows in the where χ 2 is the chi-square test value, n is the total number of samples observed, k is the lower number of categories of either variable (row or column), and V denotes the strength association between variables.
The value of V lies between 0 and 1. 57 A value near 1 indicates that the correlation is strong between variables and a value near 0 indicates a weak correlation between variables.
where d j is the difference between the two ranks of each observation,

| Relationships between polygon network types and compiled characteristics
The χ 2 tests between all geometric parameters of the TCCPs, their morphologies, and Quaternary sediments were statistically significant (P ≤ 0.05, P ≤ 0.01) and were supported with a moderate or high association between variables (≥ 0.20) ( Table 5). The intersection angles of the polygons vary according to the types of surficial   T A B L E 5 Chi-square (χ 2 ) and Cramér's V (V) values. X = the table contains cells with more than 20% with theoretical values less than five (n > 5). *P ≤ 0.05. **P ≤ 0.01. All tests were significant at a level of significance of P ≤ 0.01 and, consequently, at a level of P ≤ 0.05 sediments testify to the processes of emergence, permafrost aggradation and TCCP development (Supporting Information Figure S1).
Coastal permafrost and TCCPs are significantly younger than their counterparts found in central Nunavik, with some probably still forming on beaches and in tidal zones. 64 In the southern part of Nunavik, some TCCPs are now covered with shrubs and trees. Snow insulation associated with dense vegetation in the subarctic zone results in buffered thermal conditions that are currently not favorable for frost cracking. Drumlins, where TCCPs are mostly found, are composed of coarse-grained sediments that require colder ground temperatures for cracking (below À6 C) 65 with our observations. We also infer that permafrost adjacent to wedges located in heterogeneous coarse-grained sediments (i.e., on drumlins, eskers, and till mounds) is associated with low ground ice content due to a lower expansion coefficient and thermal stress. 67 In contrast, permafrost adjacent to smaller wedges found in fine-grained homogeneous deposits may be related to higher ground ice contents ( Figures 10 and 11).

| Morphology of TCCPs
Despite a comprehensive characterization of TCCPs over a large area, uncertainties remain regarding the polygon fills (ice, sand, or composites). However, major types of TCCPs were classified by morphology and grouped to distinguish between ice-wedge polygons, sand-wedge polygons or mixed-wedge polygons. In fact it was possible to identify low-centered TCCPs due to the trough shoulders that develop as the growing ice wedge displaces the adjacent soil. 9 Similarly, most highcentered polygons and baydjarakhs have been associated with ice wedges because they probably result from the deepening of troughs due to the melting of wedge ice along their sides (Supporting Information Figure S7). 8,9 The abundance of TCCPs in peatlands near Ungava Our data and analyses indicate that intersection angles formed according to the polygon size, local topography, and slope of the surficial deposit in which they are hosted (Table 4, Figures 5 and 6). It has been suggested that grain size distribution of soil materials and surface wetness may play a function in the development of polygonal networks and that thermal contraction coefficients may vary with different sediment deposits. 31,68,69 Thus, the stress field within certain substrates would not be homogeneous, such as in tills which are composed of various sediment sizes. 69 However, it can be deduced from our statistics that the shape and extent of the ground surface area exposed to cooling largely determine the 3D homogeneity (mainly orthogonal intersections) or heterogeneity (mainly hexagonal and random intersections) of the temperature gradients in the ground that induce the stress field responsible for thermal contraction cracking.
Our mapping demonstrates the influence of topography on polygonal network shape, where raised beaches, river terraces, flat and wide eskers, and flat and large drumlins preferentially support orthogonal networks, whereas well-drained reliefs, domed drumlins, round till mounds, and narrow eskers mainly support random and hexagonal patterns ( Figure 9).
Topography, vegetation type, and abundance of snow influence ground temperature regimes that impact polygon network geometries. local thermal stress, based on local soil moisture, snow cover depths, and vegetation conditions. 67 An example of the relationship between local topography and polygon geometries is shown in Figure 12, where surface conditions, local topography, vegetation, and snow cover were highly variable ( Figure 5).
In summary, polygon geometry is influenced strongly by local landscape topography, and relationships between slope, surface geol- For the orthogonal networks, some vertices may curve slightly over time, but the junction angles will remain mostly at 90 , and no evidence of pre-existing orthogonal cracks was found in the hexagonal polygons (Figure 5c,  This study provides the first broad-scale empirical dataset describing the spatial distribution of TCCPs, ice-wedge polygons, and their characteristics over a large part of Nunavik. The present analysis offers a better understanding of the relationship between frost cracking, surficial geology, vegetation, snow cover, and topography in Nunavik. Some TCCPs may be associated with soil and sand wedges, but additional field validation is required to further examine their respective presence and spatial distribution.
This study provides helpful data for calibrating and validating models that predict ground ice distribution in northeastern Canada. 2 The TCCP dataset will be available on NordicanaD and PINGO.

| TERMINOLOGICAL NOTE
In this paper, we refer to thermal contraction crack polygon (TCCP) as patterned ground or tundra polygons resulting from thermal stresses in the ground. It was not possible to determine with certainty from aerial photographs whether polygon sides were filled with ice, sand, or a mixture of sediments. It was more prudent to use a generic term to create a classification that includes all types of polygonal networks.
The term baydjarakh can be debated. Baydjarakh is an advanced stage of degradation of ice wedge polygon networks resulting from deepening of the troughs by melting ice wedges or by thermo-erosion. They are very similar to high-centered polygons, but baydjarakhs are more degraded, and, in some cases, the ice wedges are completely melted.
In 1988, the Glossary of Permafrost and Related terms 33 recommended against using this term and insisted on using the term "thermokarst mound." However, the term thermokarst mound removes any connotation or references to ice-wedge polygons. A comparison of baydjarakh and high-centered polygons is provided in the Supporting Information (see Figure S8).

ACKNOWLEDGEMENTS
We thank A. Boisson, S. Gagnon, and the Ministère des Forêts, de la

DATA AVAILABILITY STATEMENT
The chi-square distribution table is available in the Supporting Information ( Figure S9). Extra topography and slope models are available in Figures S10-S12. Large size maps from Figures 5 and 6 are presented in Figures S13-S17. Additional images of TCCPs are also available in the Supporting Information. The TCCP dataset will be soon available on NordicanaD and PINGO.