Advances in understanding the cooling rates and bending of needle ice: Photogrammetric and thermal observations leading to the spatial distribution of needle ice creep

The freeze–thaw cycles are a process that dominate entire landscapes with different intensity and magnitudes that are strongly changing in the current climate change conditions. In particular, the diurnal frost heave (FH) and needle ice creep (NIC) are two poorly studied processes, especially in the field, despite their effects on the sediment transportation and vegetation colonization on Alpine slopes. In a site of the Central Alps, we applied photogrammetry, thermal imaging, manual measurements and thermistors logging to improve the knowledge of the driving parameter for the formation and accretion of needle ice, of the relations of their spatial distribution with cooling rates and the frost creep modelling change in the presence of the bending process. Among the driving factors, we demonstrate that very local conditions govern the development of the needle ice, like the air cooling rates, the deepening of the freezing front and the surficial cooling rates. Especially low cooling rates are necessary for an intense process. Also the grain size influences the development of the needle ice with higher needle ice development under fine sediment than under small clasts due to a differential thermal conductivity. Moreover, the bending of needle ice has proved to be dependent on its length, age and type of heaved sediment. Therefore, potential frost creep formulas needed to be ameliorated considering the bending angle. Subsequently, a question about the role of slope gradient in NIC is treated.

In situ the range of NIC rate is very wide: It can be very low as for the Japanese Alps (0.3 cm year À1 , [Nakaya, 1995]) and very high (up to 100 cm year À1 in Bolivia [Francou & Bertran, 1997]), although the majority of the works reported values around 10 cm year À1 (Matsuoka, 1998a(Matsuoka, , 2014;;Pérez, 1987).
Close-range remote sensing techniques have been used to spatialize the effects of NIC, such as time-lapse imaging (Matsuoka, 2014;Nel et al., 2020) or photogrammetry (Ponti, Cannone, & Guglielmin, 2018).We believe that, among the remote sensing techniques, thermal imaging is fundamental to understand periglacial processes as also demonstrated at a landform scale (Ponti, Pezza, & Guglielmin, 2021).This is a key point that joins and extrapolate the processes to their climatic controls.
Numerous authors proposed formulas and corrections based on the frost heave (FH) intensity and the slope gradient (Higashi & Corte, 1971;Li, Matsuoka, & Niu, 2018;Mackay & Mathews, 1974;Matsuoka, 1998bMatsuoka, , 1998a;;Williams & Smith, 1989) to extrapolate this kind of diurnal frost creep through potential frost creep (PFC) model.However, few of them proposed alternatives to the process of the needle ice to decay (Higashi & Corte, 1971;Matsuoka, 1998a) that can largely alter the calculated frost creep.
Therefore, in this study, we want to address the following aims:

| Study area
The study area is located within the Stelvio National Park, in the Italian central Alps (Figure 1) at ca. 2670 m a.s.l.near the Stelvio Pass (2758 m a.s.l.).Here, lithology is characterized by Holocene glacial deposits and some scattered bedrock outcrops mainly siliceous (consisting of granitic and granodioritic orthogneiss, schists and biotite or two-mica paragneiss).Several periglacial landforms like solifluction lobes, terracettes, ploughing blocks and one active rock glacier (Cannone et al., 2003) are present.In the last 20 years, thermokarst phenomena were developing (Guglielmin et al., 2021) due to the permafrost degradation (Guglielmin et al., 2021).The climate is characterized by a mean annual air temperature (MAAT) around 0.0 C (closest weather station, period 2014-2020).
Precipitation values are highly variable due to the rough orography, but the total annual precipitation is generally less than 900 mm.Snow may fall at any time of the year, covering the ground continuously for more than 6 months, from November to June (Ponti, Cannone, & Guglielmin, 2021).
Permafrost has been modelled as 'mostly in cold conditions' in this area (Boeckli et al., 2012) and effectively had a patchy discontinuous distribution (Guglielmin, Cannone, & Dramis, 2001), but its lower boundary reached ca.2650 m a.s.l. on the northern slopes, while at 3000 m a.s.l.can exceed 200 m of thickness (Guglielmin et al., 2018).

| Experimental setup: Thermistors and bending
The selection of 4 plots representative of the study area includes two main different ground aspects: from 315 to 50 N for A, C and D, whereas 240 N for B. The slope was averagely 23 for A, C and D, whereas 15 for B (Ponti, Cannone, & Guglielmin, 2018).Indeed, the same plots were chosen but reduced to 50 Â 50 cm quadrats to facilitate the thermographic acquisition.The study period started on the 26 September 2018 and ended on the 18 October 2018 for a total of 23 days.This period was chosen according to Ponti, Cannone, & Guglielmin (2018) because it is the most representative of the needle ice process in this area.Close to plot C, an automatic weather station (AWS) (Onset Hobo Microstation) was installed to record at every minute the air temperature, air humidity and the soil moisture between 0 and 30 cm of depth.At the quadrat location, 3 Onset Hobo U23s (2-channel-loggers with 0.2 C of accuracy) recording every minute were placed and fixed next to each site to have 2 vertical profile temperatures at the opposite sides of each quadrat.The temperature levels corresponded to +2 cm above the initial soil surface, À1 and À2 cm below the initial soil surface.The vertical shift of the subsurface sensors was restricted using a fixed wooden slab from which only the tips of the sensor were extruded into the soil.This is in line with what previous studies suggested (Grab, 2001;Li, Matsuoka, & Niu, 2018;Nel & Boelhouwers, 2014;Yamagishi & Matsuoka, 2015).
Moreover, 5 temporally randomized soil samples (between 09:00 and 13:00 local time) were collected at 0-2 cm of depth at each plot to assess and compare the gravimetric water content (WC) variability among the quadrats.
Each ice needle growing event (INGE) was monitored prior to the formation and at the maximum growth to define the beginning of formation, the maximum accretion and the end of melting of needle ice, for each day for a total of 11 INGEs.The remaining were characterized by the absence of needle ice formation due to the snow fall (1-5 October 2018 and 15 October 2018), or due to logistic impediments, it was not possible to stay on the field at night (6-10 October 2018).
At the maximum accretion, the minimum and maximum length of ice needles (mm) at each plot was recorded with a caliper with a millimetric resolution.It is important to specify that the maximum accretion always occurs before the sunrise, while the beginning of formation can vary depending on the rigidity of the air temperature but always when the incident solar radiation at the ground surface is reduced (after the sunset).The end of melting is defined as the latest availability of solar radiation income and not as the total melting of needle ice that could happen largely before sunset.Therefore, the end of melting was always monitored at or after the sunset (ca.18:30 local time) and before the initial formation (Ponti, Cannone, & Guglielmin, 2018).
Measurement of needle ice bending has been conducted directly on the field with the same caliper.At each maximum accretion of INGE, 3 replicates of needle ice measurement were done per plot site (i.e., surrounding area adjacent to the quadrat) for a total of 113 ice needles during the study period.It measured (a) the number of growing event (cycles), (b) the type of heaved sediment (clast or fine material), (c) the length of the needles (mm), (d) the width of the needles aggregate (mm), (e) the height at which the inclusions of sediment were present in the needles (S1-5, in mm from the bottom), (f) the heights (level) at which the bending developed (B1-3, in mm from the bottom), (g) the angle between the heaved sediment and the bending level ( ) and (h) the diameter and thickness of the heaved sediment (mm) (Figure 2).

| Photogrammetry and thermography
Each quadrat was delimited with 4 stable metal rods at the vertices, and at the internal area, 5 stable toothpicks were inserted in the ground for 10-cm depth to prevent displacements.The above-ground (7 cm) top part of each toothpick was utilized as ground control point as well as for the top part of the metal rods.Because the accuracy of the best GPS would not have produced good results at such local scale, we set one metal rod as x = 0, y = 0 and z = 0 cm, and we collected 3D coordinates relative to that point 'zero' and distances (cm) among rods and toothpicks to rectify the study surface similarly to Ponti, Cannone, & Guglielmin (2018).Subsequently, 19 hemispherical and convergent (James & Robson, 2014) photographs were taken above each quadrat with a 100% overlapping pattern as described in Ponti, Cannone, & Guglielmin (2018).This repetition of images acquisition was conducted twice per INGE (11 events) that was at the end of melting and at the maximum accretion to assess the spatial distribution of the needle ice maximum length.The photographs were taken with a smartphone camera (20 MP) already used for other photogrammetric studies (Ponti, Cannone, & Guglielmin, 2018;Ponti & Guglielmin, 2021;Tavani et al., 2020) and coupled with the internal flash light during the night acquisitions.
The thermographic acquisition consisted of a portable thermal camera (FLIR E85, 384 Â 288 pixels, 0.1 C of resolution, 2.0 C accuracy) used for taking thermal images as much as perpendicularly to the plot surface from the same fixed position per each quadrat.The emissivity of this soil surface (0.97) was measured by adjusting the emissivity of the camera when the surface instant reading of a portable thermometer (0.1 C of accuracy) matched the temperature showed by the camera (Schnepfleitner et al., 2016).The atmospheric radiation was measured using the crumpled aluminum foil method described by Wolfe and Zissis (Wolfe & Zissis, 1978).
The aim was to detect the spatial distribution of the surface temperature at each development stage of INGE (end of melting, beginning of formation and maximum accretion).Only at 6 of 11 INGEs the 3 stages were visually distinguishable on the field, in particular when the beginning of formation was sufficiently slow to be individuated.Therefore, all the photogrammetric and thermographic data were considered for the 6 complete INGEs.

| Data processing
To assess the inter-plot variability of the WC during the study period, the daily minima, maxima and standard deviation were calculated among plots.The thermistors placed at each quadrat were hourly averaged and the mean for the 2 replicates calculated in order to have a unique plot-scale hourly temperature at +2, À1 and À2 cm.Subsequently, the freezing front depth (FFD, mm) at the ice nucleation was extrapolated with the ground temperature lapse at À1 and À2 cm.
The ice nucleation usually happens when a radiative cooling is present at the surface (Branson, 1992), and therefore, it should happen between the sunset and the initial formation.Here, we reasonably assumed the ice nucleation to be 1 h before the beginning of formation.In this way, it was possible to calculate the warming (+) or cooling (À) rates (WCRs, C h À1 ) at each level (i.e., air +2 cm, surface [via thermography] and À1 cm) considering the INGE stages: WCR1_S (surface), WCR1_A (air at +2 cm) and WCR1_D (at À1 cm) have been defined for the temperature change between the end of melting and the ice nucleation, whereas WCR2_S (surface), WCR2_A (air at +2 cm) and WCR2_D (at À1 cm) between the initial formation and the maximum accretion.Both a warming and cooling rate were recorded, probably due to the water changing phase from liquid to solid that slightly altered the surface emissivity of the materials in sub-plot areas (Zhang et al., 2010).WCR also depended on the position of the FF, and needle ice formation/accretion was possible because at few mm of vertical distance, both cooling and warming rates can exist.
The photogrammetric process was operated in Agisoft Metashape 1.5 and aimed to build digital elevation models (DEMs) of the plot surface with reduced errors deriving from the bundle adjustment as suggested by guidelines (James et al., 2019).To minimize the errors, the same 3D relative coordinates (cm) of the 9 ground control points (GCPs) were used for each survey, and the camera distortion was estimated (self-calibration) at the first processed model, and a precalibration of camera parameters was applied to the successive models.Indeed, the pre-calibration of cameras in the bundle adjustment diminishes the doming-effect (Griffiths & Burningham, 2019;James, Robson, & Smith, 2017).This method yielded a minimized root mean square error (RMSE) of tie points (px), GCPs (px) and distances (cm).The generated DEMs provided 1 mm of horizontal resolution, and their validation was possible by comparing the length of ice needles (minima and maxima) recorded on the field.The first days generated orthophotos were also considered to understand the distribution pattern of the clasts and fine material surface domain.
The thermal images were converted into raster format through the software FLIR Tools ® and imported into ArcGIS 10.8.Here, the thermal rasters were georeferenced with the 'projective transformation' method onto the DEMs by using the 4 metal rods at the vertices as GCPs, and the final resolution resulted to be around 2.6 mm.Always in ArcGIS, the DEMs and thermal rasters were processed in different ways.The formers were subtracted, that is, maximum accretion minus the previous end of melting to obtain the spatial distribution of ice needle lengths in one INGE; the latter were treated like the 3 stages used for the thermistors analysis in order to obtain the spatial distribution of surface warming/cooling rate (WCR1_S and WCR2_S).
For each plot a matrix with data concerning ice needle lengths, thermographic WCR1/2_S, thermistors WCR1/2_A and WCR1/2_D were produced.This matrix was treated in STATISTICA ® with a multiple regression analysis and backward stepwise option.
For a spatial analysis, at each quadrat, 9 circle polygons (best compromise between statistical representativeness and distribution homogeneity in the quadrat) with different diameters and well distributed were drawn in ArcGIS in order to extract the average ice needle lengths, WCR1_S and WCR2_S per single INGE.After classifying the ice needle length in 3 classes (1 = 0-0.5 mm, 2 = 0.5-1 The cumulative FH and the mean WCRs_S for the entire period were also divided into 2 classes (at the median value) to conduct an areal spatial distribution of the processes according to the surficial domain, expressed as clast domain or fine material domain.
The data matrix of bending observations was treated with a oneway ANOVA to understand the variability of B1 level and the bending angle according to the other characteristics.Moreover, a linear regression and RMSEs were chosen to test the accuracy of the modelled NIC implemented with the bending effect with measured displacements, whereupon compared with other theoretical creep models (Higashi & Corte, 1971;Mackay & Mathews, 1974).The observations relied on the measuring of the displacement of 15 visible clasts of different sizes in plot A for the INGE of 17/10/18 (one of the most visible events) in ArcGIS.

| Needle ice development
Generally, the maximum accretion of needle ice occurred just before the sunrise, and depending on the availability of direct solar radiation due to the clouds, it happened between 07:00 and 08:00 local time.
The beginning of needle ice formation commenced between 22:40 and 23:30 local time except on 16/10/18 when it occurred earlier (20:40 local time).This was due to the lower daily mean air temperature on that day (1.17C) compared with the rest of INGE days (2.4-8.7 C), data not shown.The end of melting was recorded between the sunset and the beginning of formation and generally laid between half an hour (29/09/18) and 2 h and 10 min earlier than the beginning of formation (for instance on 26/09/18).In this period, the recorded INGEs accounted for 6 events, and a summary of the development stages is shown in Table 1.
The topography and texture (having a mean gravel percentage of 60% for A, B and C, whereas 81% for D; and an average fine material percentage of 37% for A, B and C, whereas 15% for D) of the four plots remained unvaried from 2016 to 2018, except of the vegetation cover that was avoided in this study by selecting the vegetation free quadrant of the original 1 m 2 plot.Thus, we suggest seeing the plots' locations characteristics in Ponti, Cannone, & Guglielmin (2018).Apart from the soil temperature, the other inconstant environmental parameter was the soil WC.Despite the snowfall events and the different aspect of the plots' sites, the soil WC at 0-2 cm of depth did not show a great variability among the sites.Indeed, the standard variation of soil WC ranged between 1.2% and 6.7%, namely, indicating a WC between a minimum of 12.0% and a maximum of 29.7% that is in line with what found in Ponti, Cannone, & Guglielmin (2018) for needle ice development in the same area.Indeed, the variability of needle ice development in this study was not subjected to a large difference in soil WC (Table 2).
Even though the variability of the WC was not high, it is still fundamental for the cryosuction occurrence.Therefore, the relation between the upper soil WC and maximum needle ice length should be remarkable.However, as shown in Figure 3, neither within nor nearby the plot the amount soil moisture affected the development of needle ice in this study (R 2 < 0.1).
Because tall ice needles rarely formed during the study period (the most remarkable event occurred in A on 14/10/18 with an average length of 32.6 mm), it has been preferable to show the cumulative needle ice FH as the sum of the 6 INGEs to understand the FH spatial T A B L E 1 Times of the occurrence of the 3 INGE stages.No data (À) refer to the beginning and ending of the experiment or to the moments during which the beginning of formation was not detectable.Please note that, to facilitate the reading, cells belonging to the same color refer to the same INGE cycle.distribution (Figure 4).However, single INGE FHs were considered for the relations between daily FH and cooling rates.It is evident how and the FH minima occurred in B and C (0 cm).This fact is in line with the maximum spatial range of cumulative FH occurred in A (11.1 cm), followed by C (10.8 cm), B (7.2 cm) and D (7.0 cm).

| Types of WCRs
Overall, the northern aspect of plot A showed the deepest average FFD (À15.4 mm) and the tallest needle ice (12.4 mm).The southern aspect of plot B, instead, showed the greatest average cooling rate at FFD (À0.36 C h À1 ).
Averagely, the FFDs were shallower or even above the ground level at the beginning of needle ice formation, whereas deeper during the maximum accretion, around at 1 cm of depth.Therefore, it is interesting to show the FFDs and the WCR at 1 cm of depth.The greatest daily FH occurred in A (25.2 mm), followed by C (18.9 mm), D (17.0 mm) and B (16.5 mm), even though the maximum FFDs did not follow the same pattern: A (À49 mm) followed by B (À44 mm), C (À42 mm) and D (À32 mm).The lowest negative WCR2_D occurred in D (close to À0.02 C h À1 on 18/10/18) and in A (close to À0.03 C h À1 on 18/10/18), followed by C (À0.04 C h À1 on 18/10/18) and B (À0.07 C h À1 on 18/10/18).Only in A the lowest WCR2_D matched the maximum FH among the plots, and additionally, it did not follow the maximum FFDs among the plots (Table 3).A statistically significant negative linear regressions were found in A between FHs and WCR2_D (R 2 = 0.66, p < 0.05) and between FHs and FFDs (R 2 = 0.82, p < 0.05) (data not shown).4).
If we consider all the subsample polygons at all the plots for the whole period, we can extract the spatial mean values of single WCRs.
These values, ordered in 3 classes of needle ice length (1 = 0-0.5 mm, 2 = 0.5-1.0mm, 3 = >1.0mm), show different trends from small to high needle ice.Indeed, the more the WCRs got close to 0, the highest needle ice was observed.This is particularly true for WCR1_S with values ranging from À0.49 and 0.10 C h À1 and for the sum of WCR1_S and WCR2_S with values ranging from À0.55 to À0.11 C h À1 .Conversely, for WCR2_S, the trend was opposite and less evident because of the lower variability of the WCRs that ranged from À0.12 to À0.18 C h À1 .The resulting ANOVA among the groups of needle ice length showed the means statistically differed, with less robust results for WCR2_S (Figure 5).

| Spatial variability of surface WCRs
The average spatial distribution of WCR1_S and WCR2_S was different per each plot.It is important to notice that not only a cooling rate was responsible for the beginning of formation but also a warming T A B L E 2 Inter-plot variability of the water content (%) for the samples collected at 0-2 cm of depth.Similarly, before the maximum accretion, A showed the cooler WCR2 (mean of À0.19 C h À 1), followed by D (mean of À0.14 C h À1 ), then B (mean of À0.11 C h À1 ) and C (mean of À0.08 C h À1 ).In addition, generally, the WCRs were greater at the first development stage rather than before the maximum accretion (Figure 6).
By focusing on the distribution of the total heave under the fine material surface or clasts, it is observable that there is not a unique trend and, especially, that not always the clast domain showed a lower cumulative heave.The largest mean heave difference was recorded in plot A (4.1 mm) showing a smaller heave for the clast domain, whereas the smallest was recorded in plot B (0.9 mm) showing a smaller heave for the fine material domain.More interestingly, the plots that undergo larger total heave were those with a dominant fine material surface cover, suggesting that non-concentrated clasts were favored by FH (Table 5).

| Bending process
Many morphometric parameters were observed during the maximum accretion of needle ice in order to highlight the bending process.On 113 needles investigated, 82 showed at least one inclusion of sediment (S1), whereas only 9 showed a second bending level (B2).Generally, the first inclusion of sediment level was lower (8 mm) than the bending level (14 mm) as well as for second or third inclusions and bending levels.The median of needle ice length was 26 mm and the width of needles' aggregates relatively small (7 mm), but this is dependent on the diameter of the topping sediment that induced a cohesion of needle ice at its bottom.The diameter of the topping sediment type (clast or fine material) was not very large (5 mm, overall) but its thickness greater (14 mm).More importantly, the median of bending was 20 with a maximum value of 90 that is remarkable for the frost creep process (Table 6).An interesting result was the minimum limit for needle ice to bend: It was observed to occur at 2 mm from the needles' bottom.
Because the first level of the bending (B1) and the bending angle are the most important parameters to understand the NIC, two box and whiskers plots are presented to remark B1 and angle differences among the morphometric parameters separated by their median value.
Concerning B1, different distribution patterns are visible for S1, bending angle and length.In contrast, concerning bending angle, visible differences occur for B1 and the type of topping sediment (Figure 7).
T A B L E 4 Summary of the statistical tests of each dependent variable that affected the needle ice length in the multiple regression analysis.
Please note that symbols refer to direct (+) or inverse (À) proportionality between the investigated variable and the needle ice length.When looking at the statistical analysis of the distribution patterns of the morphometric parameters (Table 7), it is clear that S1, bending angle and length are statistically well discriminated by B1.
Similarly, B1, length and sediment type (clast or fine sediment) are statistically well discriminated by the bending angle.In other words, the length of needle ice statistically affects both B1 and the bending angle, whereas S1 affects B1 and the sediment type the bending angle (Table 7).

| Needle ice development
Our findings suggest that there is an inter-plot variability that is crucial for determining the length of needle ice and in turn the subsequent creep magnitude.Indeed, at sites A and C, the cumulative FH resulted to be higher than at the sites B and D. As a confirmation, we found the deepest average FFD at A (15.4 mm) that was northexposed, whereas the shallower FFD at D (7.7 mm) that was W-exposed.B and C showed shallow FFD, as well, with an average of 9.9 and 9.1 mm, respectively.Moreover, it was shown that the upper soil WC does not affect the maximum length of needle ice, probably because the cryosuction could be more or less effective independently from the soil moisture.
The reduced surface energy balance at northern aspect provides more INGEs because on 18/10/18, plot B (and D) did not show a detectable growth of needle ice.This is in line with Pérez (1987) but the opposite to what found by Ponti, Cannone, & Guglielmin (2018) despite the same locations.It is true that aspect and thus the radiative balance could determine the diurnal FH (Lawler, 1993) T A B L E 6 Summary of bending measurements and the other morphometric parameters described in the methods.S1-5 refer to the distance from the ground surface to the sediment inclusion, and B1-3 refer to the distance from the ground surface to the bending occurrence.temperatures.Rather, as well demonstrated by other authors, it is a matter of cooling rates (Grab, 2001;Lawler, 1993;Outcalt, 1971;Ponti, Cannone, & Guglielmin, 2018;Yamagishi & Matsuoka, 2015) that are controlled not only by radiation.

N of cases
In support of this, the only WCR that weakly governed the length of needle ice through the plots was WCR1_A, similarly to what found by Ponti, Cannone, & Guglielmin (2018) in their model.Surface or ground cooling rates, instead, were key drivers at the plot scale that depended on local topo-edaphic characteristics.Therefore, the plots underwent a differential process probably more linked to the possibility of cryosuction than to the weather conditions.Therefore, the surface air temperature cooling could be considered as a common trigger of needle ice development, and ground level-specific WCRs and/or thermal conductivity of the surface provided a differential cryosuction for the FH occurrence (Meentemeyer & Zippin, 1981;Outcalt, 1971).

| Warming rates
We observed few warming rates during the INGEs as also found by Grab (2001) and Nel et al. (2020).This is explained by the sensible heat transfer from deeper and warmer depths to the cooling surface after the sunset showing that the ice nucleation could happen at the very top millimetres (Ponti, Cannone, & Guglielmin, 2018).Positive values also occurred for WCR2_S, but it was related to the sensible heat transfer of warmer air temperatures that can happen in absence of nightly breeze or to the emissivity change during the water freezing in WCR1_S (Zhang et al., 2010).Indeed, the wind action usually hinders the development of needle ice at the ground surface by changing the vertical fluxes of energy (Grab, 2001) or, more likely, drying the very top ground surface.More interestingly, warming rates at 1 cm of depth were only recorded in A and C and only at the formation phase.
Indeed, positive WCRs at a negative FFD mean that an increase of ground temperature happened but still remained below 0 C, that is, a slowing down of cooling rate that indicated ice nucleation (Lawler, 1993;Outcalt, 1971).Further, positive WCRs with above ground FF are possible at the beginning of formation because the FF is still approaching the surface before the visualization of the first ice needles.WCR2_A has been the driver parameter common to each plot (and independent from the ground aspect) that showed a good relationship with needle ice accretion except in C.This is related to the fact that the sensible heat transfer from the ground surface to the air is not as much remarked as through solid means.Therefore, despite the air temperature measurements were close to the ground surface, the air cooling rates reflected the stability of the ground boundary layer (Barry & Chorley, 2021) that has been similar among the plots.Moreover, this stability occurred at the accretion phase, so when air temperatures slowly dropped with low cooling rates.Indeed, many authors found that high cooling rates prevent the development of needle ice once formed (Grab, 2001;Lawler, 1993;Outcalt, 1971;Ponti, Cannone, & Guglielmin, 2018;Yamagishi & Matsuoka, 2015).In this study, we confirm the development of needle ice under low ground cooling rates, well below than the ones found by Grab (2001) (2.4 C h À1 ).However, not always lower cooling rates corresponded to taller needle ice.

| Spatial variability of surface WCRs
Our intra-plot and multi-temporal data demonstrate that the surface spatial variability of WCRs and needle ice length match.Indeed, relatively high (À0.49to 0.10 C h À1 ) surface WCRs are needed to form needle ice (Lawler, 1993;Nel & Boelhouwers, 2014) but not extremely high, otherwise ice lenses form (Branson, 1992;Higashi, 1958).On the other hand, lower surface WCRs (À0.12 to À0.18 C h À1 ) maintain a continuous accretion once formed (visible at the surface) (Figure 5).This is interesting because, as occurred at the average plot level, higher cooling rates firstly permit the FFD to deepen from the top surface and afterwards low cooling rates (stable FFD) after the formation can provide a continuity of water suction and thus ice to segregate (Lawler, 1993;Outcalt, 1971).
What joins plot B and D is the more homogenized surface cover percentage of fine material (50.9% and 52.1%, respectively) and clasts instead of A and C (84.9% and 69.4% of fine material, respectively).
This distribution provided a more homogenized surficial thermal conductivity that resulted in similar WCR_Ds (A to C and B to D) and similar total heaves per surface domain.Indeed, it was discussed that the higher thermal conductivity of stones compared with soil can provide an earlier ice nucleation (Branson, 1992) or faster FF deepening and thus little accretion.On the contrary, A and C showed a generally slight increase of total heave under the fine material domain compared with clasts and also averagely higher than B and D. This is in accordance to Li et al. (2021) and Yamagishi & Matsuoka (2015) that demonstrated that taller needles develop in fine material and move clasts towards the clast domain, increasing their concentration.We additionally claim that a well-mixed surface cover hinders the spatial autocorrelation (homogeneous patches) of thermal conductivity and thus the formation of tall needle ice.

| Bending process and frost creep
Our findings demonstrated that the bending of needle ice could happen at levels such as 2 mm from the ground and with a median of 20 .This is the first time that in situ limits of bending have been recorded, T A B L E 7 Summary of the one-way ANOVA conducted on B1 (mm) and bending angle ( ), grouped by other morphometric parameters.Each grouping parameter is divided in 2 classes: lower and greater than the median value (showed in parentheses) or according to the aspect of the ground or the type of topping sediment.similarly to what was defined for toppling (Nel et al., 2020;Ponti, Cannone, & Guglielmin, 2018).This is important because the effects on frost creep might be affected by almost the total length of needle ice as well as the bending angle.
In particular, what defines higher B1 and bending angle is the length of single ice needles that suggest a greater 'fragility' probably due to their rapid growth or more growing events (polycyclic) compared with smaller ones.Indeed, the ANOVA showed that taller needle ice has a higher number of single-night cyclic growths ( p < 0.005) (data not shown) that indicate a hydric stress during formation (Lawler, 1993;Soons & Greenland, 1970).Further, the greatest bending angles were found at higher B1 levels probably because additional bending was possible at greater distance from the ice nucleation level, so at the 'older' and more fragile part of the single ice needles.
The inclusion of sediment during the needle ice growth was already observed and also treated as an indicator of rapid accretion when the FFD descends during higher cooling rates (Branson, Lawler, & Glen, 1996;Lawler, 1993).In this study, S1 was generally lower than B1.This demonstrates that bending might occur: (a) with relatively low cooling rates during the growth (as opposed to the inclusion of sediment), (b) generally (76%) earlier (higher B1) than the last inclusion of sediment (S1) and rarely (20%) at the same time (level).This consideration suggests that the inclusion of sediment is not responsible of the bending, neither through the different thermal conductivity nor through the acceleration of cooling rates.Therefore, the fact that higher S1 corresponds to higher B1 could only be related to the length of needle ice so its growing rate.Lawler (1993) associated the bending process to the wind action.
However, in our study, we only met a gentle night breeze and a random direction of bending that made us exclude the wind hypothesis.
Rather, we found that ice needles with topping clast were less bent than needles with topping fine sediment.This is likely due to the fact that topping clasts were associated to shallower ice nucleation, thus a thinner layer of fine material between the clast and the nucleation point.In support of this, the absence of fine sediment between the heaved clast and the ice needle suggests a very surficial development compared to the fine sediment heave.
As a matter of fact, we found that the maximum thickness of where h is the length of ice needles, when the creep magnitude is mainly given by the toppling process (Mackay & Mathews, 1974).
When the creep magnitude is attributed to both toppling and bending on a slope gradient (Higashi & Corte, 1971) as where h is the length of ice needles, θ the slope gradient and C a constant coefficient of particle cohesion that Matsuoka (1998b) found to be equal to 3 to better explain NIC from FH.However, in this study, we measured the bending process in situ, and thus, we suggest an alternative to the previous equations based on the measured bending angles.Because toppling is more relevant than the slope angle from field observation (Mackay & Mathews, 1974), the bending contribution to the NIC (C b ) is thus given by where θ is the slope angle, h is the length of ice needles, B1 is the height of the bending level and α the bending angle (Figure 8).
At a first glance, there is no relation between the observed and modelled (Equation 3) creep due to the very low R 2 value.However, this relation does not consider the direction of the creep that, in turn, affect the intensity.Indeed, purely based on the NIC magnitude, the comparison of the 3 models shows that the model developed according to Equation 3 has the least RMSE (4.4 mm) with the observed displacements, followed by Equation 2(5.1 mm) and then Equation 1 (15.5 mm) (Figure 9).The PFC formula that considers all the frost creep processes (Matsuoka, 2001;Washburn, 1979) was used in literature and many corrections for NIC proposed (Higashi & Corte, 1971;Li, Matsuoka, & Niu, 2018;Mackay & Mathews, 1974;Matsuoka, 1998a).However, it still remains difficult to forecast the NIC under the effect of the bending process.Others found it to be dependent on the type and the spatial cover of needle ice (Branson, 1992), and lab experiments demonstrated how bending mainly occurs at the soil/clasts domain border (Yamagishi & Matsuoka, 2015).
In this study, we met the necessity to measure the bending process that has never been monitored on the field before.Hereby, the proposed formula divides the final creep in 2 components: the straight part of the ice needle (B1 level) and the contribution of the bent part (length minus B1) (Figure 8).On the basis of our field investigation, we sustain that Mackay & Mathews' (1974)  more than the slope gradient as we always found randomly bent needles on slopes of (15-24 ) (Figure 10), oppositely to Li, Matsuoka, & Niu (2018).Indeed, because the tendency of decay is the toppling process because melting also occurs at the base of needle ice (Mackay & Mathews, 1974;Matsuoka, 1998b;Yamagishi & Matsuoka, 2015), the gravity force induces the direction of decay of the randomly bent ice needle with a topping sediment not considering the slope direction.Therefore, we do not agree with Li, Matsuoka, & Niu (2018) and Matsuoka (1998b) because both on gentler or steeper slopes, the real NIC is determined only by the direction of bending and not by the slope gradient.A rare exception would be a combination of high uphill bending angles (max 90 ) on such very steep slopes (i.e., river banks) that the gravity force prevailed downhill in any case.

| CONCLUSION
In this study, we explained the role of different types of cooling rates in needle ice development.Moreover, we recorded data concerning the bending process in situ for the first time and proposed a correction of the PFC models that rarely considered the variability in the field.In fact, very local conditions (WC, snow, thermal conductivity, cooling rates) drive the development of taller or smaller needle ice.
Generally, air cooling rates favor the ice nucleation, even though ground warming can occur during the formation and accretion due to sensible heat transfer from warmer ground layers to the cooling surface.Definitely, the ground aspect and presence of soil moisture determine the shift of the FF that, with specific cooling rates, promotes the accretion of needle ice.In particular, relatively high surface cooling rates are needed for the formation of needle ice to reach deeper moisture, whereas lower surface cooling rates with a stabilized FFD yield the ice segregation providing taller needle ice.On a spatial scale, the distribution of clasts is able to affect the intensity of the FH because of the higher thermal conductivity of rocks that deepens the FF fastly, avoiding the formation of taller needle ice.Indeed, only where fine sediment is well visible, a good accretion is noticed.
Surprisingly, the bending process, that largely affects the frost creep with an average of 20 , can start at 2 mm from the bottom of needle ice.Field measurements demonstrated that the bending level (B1) and angle are a consequence of the fragility of the single ice needle when tall (owing to the rapid or polycyclic accretion) and thus having an old top portion, as well as the type of heaved sediment.
Conversely, the inclusion of sediments within the ice needle does not affect the bending angle.Topping clasts are also an indicator of shallower ice nucleation that produces little angles compared with the fine sediment.Indeed, fine sediment addresses a nucleation that reaches 17 mm of depth and more intense bending.In this study, we also proposed a new PFC formula that considers the bending process that turns into creep independently from the ground slope.Rather, the creep is triggered by the gravity force acting on the weight of the heaved material, randomly bent (directed).
These considerations arise 3 scientific questions that still need to be answered to fully understand the spatial distribution of NIC: Is it possible to forecast the direction of bending needle ice?What is the areal pattern of the bending process?Does really matter the ground slope in NIC?
(a) to understand at what level the air/soil cooling rate contributes to the needle ice development; (b) to assess the spatial variability of the surficial cooling rate as potential contributor to the NIC; (c) to deeply analyze the bending process of needle ice in order to correct the current PFC models.

F
I G U R E 1 Settings of the needle ice plots that match the names in the paper (a-d) and their location close to the Stelvio Pass area (e).(f) shows a detail of plot B with the location of the temperature loggers.[Color figure can be viewed at wileyonlinelibrary.com] .0 mm, F I G U R E 2 Closeup of needle ice with indication of some processes occurred during the development and monitored in this study.The width of the heaved clast is 20 mm.[Color figure can be viewed at wileyonlinelibrary.com] 3 = >1.0mm), an analysis of variance (ANOVA) was run on this dataset in order to find the relations between growth and WCRs.
plot A was subjected to a greater cumulative FH, followed by plot C, plot D and plot B. The spatial distribution of the FH gradient was particular in A, showing a higher magnitude in the southern central part, while similar for B, C and D, showing a lateral gradient, less evident in D. Among B-D, B and C showed the highest magnitude heave at the western side, whereas D at the eastern.Moreover, the cumulative FH maxima were recorded in A and C (14.0 and 10.8 cm, respectively), At a detailed scale, considering the intra-plot spatial distribution of needle ice length and all types of WCRs, interesting findings can be observed.Generally, the most robust statistic model of different needle ice lengths occurred in D (R 2 = 0.99), followed by A (R 2 = 0.87), C (R 2 = 0.81) and B (R 2 = 0.74).By looking at the variables that more largely affected the needles lengths (high F values), it is evident that they are different at each plot.However, a common WCR2_A is found for A, B and D. This underlines the importance of temperature changes above ground also at the accretion phase of needle ice.C implemented this effect with a large contribution of WCR1_A (F = 32.8)instead, thus indicating a major contribution of the formation phase.However, A and C were subjected to the variability of WCR_D in different phases, whereas D to the WCR1_S of the formation phase (F = 89.6)(Table

F
I G U R E 3 Linear regression between the upper soil water content (%) close to the plot and the maximum length (mm) of needle ice within the plot (red) and nearby the plot (black).[Color figure can be viewed at wileyonlinelibrary.com] rate.This is true especially for plot A, where the warming rates reached 1.68 C h À1 , whereas for B, C and D reached lower values: 0.87, 1.28 and 0.84 C h À1 , respectively.Conversely, between the beginning of formation and the maximum accretion, the maximum warming rate of A resulted lower (0.01 C h À1 ) than the others (B = 0.29, C = 0.19, D = 0.17 C h À1 ).Further, plot A was the only location at which the temperatures change displayed a high spatial autocorrelation, whereas the other plots showed a low spatial autocorrelation.Moreover, the patterns of WCR1_S and WCR2_S in A were generally opposite, showing that what cooled rapidly at the beginning of formation slowed down during the accretion and what warmed rapidly at the beginning of formation then cooled rapidly during the accretion.A similar pattern is visible in B. The spatial autocorrelation pattern of WCRs in A showed a similar distribution for its needle ice length as well as slightly in B, whereas in C and D the high spatial autocorrelation of needle ice length did not follow the low spatial autocorrelation of their WCRs (Figure 6).Averagely, during the beginning of formation, plot A showed the warmer WCR1 (mean of 0.55 C h À1 ), followed by D (mean of À0.3 C h À1 ), C (mean of À0.47 C h À1 ) and B (mean of À0.77 C h À ).

F
Abbreviations: INGE, ice needle growing event; WCR, warming or cooling rate.

F
I G U R E 5 Mean values and 95% confidence intervals of WCR1_S (A), WCR2_S (B) and the sum of the two previous WCRs (C) for the 3 classes of needle ice lengths (1 = 0-0.5 mm, 2 = 0.5-1.0mm, 3 = >1.0mm) recorded in different days within the 9 randomly distributed polygons in the quadrat.The ANOVA results are represented with the same letters of the panels.[Color figure can be viewed at wileyonlinelibrary.com]

F
I G U R E 6 Intra-plot spatial distribution of the average surficial cooling rates from the end of melting to the initial formation (WCR1_S) and from the initial formation to the maximum accretion (WCR2_S).(a) WCR1_S for A, (b) WCR2_S for A, (c) WCR1_S for B, (d) WCR2_S for B, (e) WCR1_S for C, (f) WCR2_S for C, (g) WCR1_S for D and (h) WCR2_S for D. [Color figure can be viewed at wileyonlinelibrary.com]

F
I G U R E 7 Box ad whiskers representation (non-outlier min, max, 25/75 percentiles and median) of the height of bending (B1, mm) and the bending angle ( ).Variables are grouped according to the median value or according to the ground aspect or the type of topping sediment.[Color figure can be viewed at wileyonlinelibrary.com] topping fine sediment is 17 mm that corresponds to the maximum depth of the ice nucleation level in this study, and it is shallower to the displacement depth found byMackay & Mathews (1974) (8 cm), Pérez (1988) (3-5 cm) and Nakaya (1995) (4 cm) because they also considered ice lensing.However, our thickness is similar to what found by Yamagishi & Matsuoka (2015) (1 cm), which solely focused on needle ice.Considering the existing PFC formulas, previous research defined the theoretical NIC as

F
I G U R E 8 Scheme representing the creep of a downward bent ice needle: (a) straight ice needle, (b) bent ice needle, (c) toppling of a bent ice needle, (d) zoom of the fallen ice needle with the final components (x and y) of the creep and the measured variables on field.[Color figure can be viewed at wileyonlinelibrary.com] formula, which considered the prevalence of the toppling process and excluded the slope gradient, was very far from our results.Oppositely,Higashi & Corte's (1971) formula, which considered both toppling and bending processes with the slope gradient, gets close to our creep results.However, our proposal best fits the observed values probably because the bending process is much more frequent than expected and needed to be quantified and applied.The poor statistically result between observed and modelled creep could be explained by the effect of the direction of NIC.As observed, the bending direction is random and ultimately affects the fall of the ice needle influencing the final creep indipendently from the ground slope.Therefore, an uphillbent ice needle will fall upslope and give a discrepancy between the proposed model and the observed creep.Likewise, a creep could occur in perfectly plan surfaces due to bending.Therefore, as a consequence of the negligibility of the slope angle, we do not agree with the correction of PFC formula made byLi, Matsuoka, & Niu (2018) that underestimated the NIC for few runs of the lab experiment.Indeed, bending of needle ice in field could both underestimate (downslope curvature) or overestimate (upslope curvature) the NIC depending on the decay direction.According to us, the weight of the topping material on a bent ice needle affects the direction of decay F I G U R E 9 Linear regression between observed and modelled needle ice creep (NIC) through Equation 3 for the ice needle growing event (INGE) of 17/10/18 in plot A. Please note that an average of B1 = 27 mm and α = 20 were used in Equation 3. On the right, the RMSEs between observed and modelled NIC using the different potential frost creep (PFC) formulas are shown.F I G U R E 1 0 Example of random bending (uphill and downhill) on gentle slopes adjacent to the study plots.The bag in (a) is 15 cm long, whereas needle ice in both (a) and (b) are approximately 8 cm tall.[Color figure can be viewed at wileyonlinelibrary.com]