Quantification of freeze–thaw hysteresis of unfrozen water content and electrical resistivity from time lapse measurements in the active layer and permafrost

We quantified seasonal, in‐situ variation in ground electrical resistivity, unfrozen water content, and ground temperature using daily and sub‐daily measurements in active layer and in permafrost on a high‐latitude monitoring site in West Greenland. The in‐situ unfrozen water content in a clayey soil undergoing phase change was 10% higher during freezing than during thawing at the same ground temperature, thus shifting the freezing curve parameterization seasonally. The freeze–thaw unfrozen water content hysteresis was a measurable part of the energy balance of the active layer and corresponds to 12.1% of the total yearly enthalpy change for our field site. We also observed and quantified hysteresis in the relationships between ground electrical resistivity and unfrozen water content, where the inverted resistivities during freezing in the temperature range from 0°C to approximately −4°C were up to one order of magnitude higher than during thawing at the same unfrozen water content. As a result, the petro‐physical models that uniquely derive soil unfrozen water content from bulk soil electrical resistivity measured from ground surface do not fully capture the complexity of the in‐situ processes. The hysteresis should therefore be considered in quantitative interpretation of ground resistivity changes in terms of ground temperature and ice content changes.


| INTRODUCTION
As permafrost is a temperature-defined phenomenon, ground temperatures are the most commonly cited indicator of permafrost conditions.Yet in case of warm (with mean annual ground temperature [MAGT] near the thawing point), fine-grained, and/or saline permafrost, ground physical and mechanical conditions cannot be fully ascertained from the temperature observations alone.For example, long-term permafrost temperature records show twice the rate of increase in MAGT in the colder, continuous permafrost zone in comparison with the warmer, discontinuous permafrost zone. 1 The lower rate of temperature increase in warm permafrost does not, however, imply more stable ground conditions, as gradual phase change happens at zero and sub-zero temperatures.For a number of practical applications-including permafrost thermal modeling, long-term predictions and engineering considerations-it is thus the unfrozen water content that determines the bulk soil thermal properties 2 and their mechanical behavior. 3Soils in permafrost landscapes typically spend well over a quarter of the year in the temperature range between À5 C and 0 C 4 where the bulk of the phase change processes occurs.
Knowledge of the soil-specific freeze-thaw regime is thus critical for understanding of permafrost processes, for realistic predictions of thaw degradation potential 5 and associated ground stability changes. 6ng-term monitoring is an essential tool to arrive at this understanding. 7frozen water content monitoring in-situ is typically done using time-or frequency-domain reflectometry sensors (TDR or FDR) (e.g. 4,8,9); however, extensive field data sets from high-latitude permafrost remain sparse.
Electrical resistivity tomography (ERT)-a noninvasive geophysical method highly sensitive to phase change between water and ice-has been successfully used to monitor temporal, relative changes in unfrozen water content in high-latitude permafrost, 10 in mountain permafrost, 11 and rock walls. 12Studies by 13 and 14 show that there is a quantitative link between electrical properties and temperaturedependent moisture content of geological materials.Calibration of this link requires thorough understanding of petro-physical relationships and thus availability of data with sufficient level of detail, spanning several freeze-thaw cycles.6][17] Although the bulk ground resistivity depends on a number of factors (including soil mineralogy, porosity, fraction of unfrozen pore water, and geochemical composition of the pore water 18,19 ), the aqueous solution is, in most cases, the main conducting phase in soil.The aqueous solution and the ice content are also the only phases that substantially change volumetric proportions over the course of a year.
In this work, we present 3-years-worth of data that map the relationships between ground temperature, soil unfrozen water content, and electrical resistivity at a high-latitude permafrost monitoring site in Ilulissat, West Greenland.Comparatively dense acquisitions and joint interpretation of three types of data provide an unprecedented level of insight into the quantitative relationships between these soil parameters.We also highlight challenges that are associated with interpretation of data acquired in challenging environmental conditions, at different temporal and spatial scales, and the need for sitespecific calibration of the petro-physical relationships.Ultimately, these considerations should lead to more realistic modeling assumptions, improved estimates of temperature-dependent physical and mechanical soil properties, and more useful model predictions.

| ILULISSAT SITE DESCRIPTION
The permafrost monitoring site is situated near the airport in Ilulissat (69 14' N, 51 3' W, 33m above sea level), on the Greenland mainland in the inner part of the Disko Bay (Figure 1).The mean annual air temperature between 2010 and 2019 was À3.1 C (data from 20 ).The site is located in the continuous permafrost zone. 21,22e soils at the site are classified as silty to very silty clays according to the Danish engineering geological practice, 23 which is normally used in Greenland (soil profile in Figure 2a).The soils have a very high fines content (d 2μm ¼ 30 À 50% and d 63μm ¼ 80 À 95%), with a liquid limit around 30% and plasticity index around 15%.They plot above the A-line in the Casagrande diagram and are classified as lean clays according to the Unified Soil Classification System (ASTM D2487-11).
The bedrock-Nagsugtoquidian gneisses with amphibolitic bands-is encountered at 7m depth according to the site investigations for the nearby airport in 1978 (borehole ILU78020, drilled 1m from the borehole ILU2007-01). 24The soil column is covered by a few centimeterthick vegetation layers (Figure 3b).
The sedimentary profile is a result of the marine transgression that followed the deglaciation of the area around 9600 BP. 25 Isostatic rebound has since exposed the sediments to precipitation and percolation, as well as colder atmospheric conditions following the end of the Holocene optimum. 26The marine fine-grained deposits are therefore partially leached, with an upper, leached, and ice-rich permafrost layer, and a deeper, unleached zone with low ice content (as indicated by gravimetric water contents and salinity profile, see Figure 2b-d).
The active layer thickness at the monitoring site is approximately 0.9m 27  (see Figure 2d).Ice contents decrease with depth as the salinity increases to full sea water salinity (chloride contents around 19 g/L).According to an empirical formulation by 28  The depth of zero annual amplitude is 5m in ILU2013-01 (defined as the depth of maximum annual amplitude <0.1 ).
Based on this information-and supported by a number of boreholes drilled in the area between 1979 and 2018-the geological profile of the sedimentary basin can be considered having a horizontally layered structure.In summer, this consists of four layers: the active layer (expected resistivities ρ < 500Ωm), ice-rich permafrost (ρ > 5000Ωm), saline permafrost (ρ5Ωm), and the frozen bedrock (ρ > 2000Ωm).In winter, the active layer is frozen, while little seasonal variation is expected in the properties of the layers below.

| AUTOMATED MONITORING STATION
The monitoring station in Ilulissat was established in August 2012 and consists of ground temperature logging in two boreholes, an air temperature sensor, soil moisture sensors in the active layer, and an automated ground resistivity monitoring system (for overview and detailed pictures of the in-situ instrumentation, see Figure 3a-e).In the context of this article, measurement and datapoint refer to one individual measurement or datapoint, acquisition refers to the measurements or datapoints collected during one logging period (e.g.daily acquisition of resistivity measurements), and data set refers to all measurements of the given parameter collected for the duration of the study.
Air temperature was measured at the height of 1.1m above the ground surface every 3 h by a TMC-HD thermistor placed in a radiation shield and connected to a HOBO U12-008 datalogger.
Ground temperatures were measured every 3 h at two depth resolutions.Borehole ILU2013-01, drilled in August 2013, was instrumented with 3 HOBO U12-008 loggers and 12 TMC-HD thermistors at depths 0, 0.25, 0.50, 0.75, 1.00, 1.25, 1.50, 2.00, 3.00, 4.00, 5.00, and 6.0 m.In addition, detailed soil temperatures were measured in the active layer and the top of the ice-rich permafrost with a 1.5-mlong rigid thermistor string (MRC probe) with thermistors spaced by 10cm.Data from borehole ILU2007-01 (drilled in 2007 ca. 1 m from ILU78020, 27 ) were used for comparing past and recent temperature envelopes at the site; recent recordings from the borehole were, however, not used due to severe frost heave that has affected the installation in later years.
The soil unfrozen volumetric water content θ w was measured insitu by FDR, using Stevens Hydra Probe ® sensors.Two sensors were installed at the depths of 0.3m and 0.55m, respectively, and they have been logging dielectric permittivity of the soil in the active layer every 3 h.The unfrozen volumetric water content θ w ½m 3 water =m 3 bulk was calculated from the permittivities using a conversion equation suggested by Seyfried et al. 29 : where E r is the real part of the measured permittivity and the two  The ERT instrumentation consisted of an ABEM SAS1000 resistivity meter and an external electrode selector ABEM ES10-64.The ERT array was equipped with 64 electrodes, equidistantly spaced by 0.5m and individually connected to the electrode selector via coated copper cables (Figure 3c and d).Both electrodes and cables were buried at a depth of 10-20 cm to protect them from damage and removal by people and animals, and to ensure that they were installed in the clay soil below the organic surface layer.In the first year of station operation, the electrodes used were stainless steel rods (length = 8.0cm, diameter = 1.0cm, surface area = 27.0cm 2 ).In August 2013, the rods were replaced by square stainless steel meshes (side length 10.0cm, thickness 0.6cm) which, due to their increased surface area (985cm 2 ), effectively reduced the extremely high grounding resistances that hampered data acquisition throughout the winter 2012/2013. 30According to derivations of Ingeman-Nielsen & Tomaškovičová (2015, pers.comm.), the geometrical effect on the apparent resistivity measurements resulting from the finite electrodes of the size of the mesh electrodes and spaced by 0.5m is less than 1%.
In February 2014, the station was equipped with an on-site computer that has allowed for remote access and re-programming of the resistivity system, as well as daily synchronization of the collected resistivity, temperature, and soil moisture data.The data are transmitted over a cellular modem to permanent storage at a server.
The station was connected to a 230 V power outlet at a nearby Ilulissat airport facility (Figure 3b).Power was supplied to the on-site computer, the resistivity instrumentation, and the Campbell datalogger by a custom-built power supply unit that controlled the alternate supply and charging of two 12 V, 60 Ah lead-acid batteries.The batteries ensured uninterrupted operation even during power outages.The power consumption of the entire station-with the resistivity system being the main consumer-was on average 30 kWh per month in winter and 25 kWh per month in summer.These values match closely our initial power consumption estimates based on the ABEM SAS1000 specifications, which predicted the need for supplying 29 kWh of electricity in frozen, winter conditions to the resistivity system alone.
All the instrumentation was housed in a lockable aluminum box placed on a wooden pallet, to avoid disturbance to electric current paths in the ground.
The ERT system has been collecting ground resistivity measurements once every 24 h since August 2012.In this article, we focus on The filtering removed up to 8% of the Wenner-Schlumberger measurements in the first year of the station operation-while the array was equipped with the rod-shaped electrodes.After the electrode replacement by the mesh electrodes and improvement of the electrode grounding, only 0.6% of the measurements had to be discarded.

| INVERSION OF THE RESISTIVITY DATA
The Wenner-Schlumberger apparent resistivity measurements at each of the 33 vertical focus depths were averaged, with the standard deviation of the four to nine averaged measurements used as an estimate of the measurement error.This resulted in daily resistivity sounding acquisitions with up to 33 measurements each-one apparent resistivity value per each vertical focus depth.The 1D acquisitions were inverted using the AarhusInv inversion software. 32 The resistivities between 2 and 3 m (resistivity model layer 9, 1,000 Ωm) were constrained to approx.±20% (STDF = 1.2) to allow some resistivity variation due to variation in temperature at this depth.
F I G U R E 4 Initial model for the 1D inversions, consisting of 12 layers of fixed depths and thicknesses, with the 12 th layer extending to infinite depth (lower half-space).The choice of layer boundaries was based on geological information from the boreholes ILU78020 and ILU2007-01 (see Figure 2).The black dots indicate the vertical focus depths (depths of maximum sensitivity) of the Wenner-Schlumberger quadrupole configurations.
The a priori constraints applied to the layers 1-9 (Figure 4) were not fixed limits, but rather included in such a way that if the data strongly favored a larger difference, the model could accommodate it. 33) Use of the loose vertical constraints reduced overfitting of the data, as justified by the sensitivity analysis (5.7).
To fix the resistivities in layers below 3m to plausible values, Schlumberger soundings from the same location from years 1979 34 and 2007 35 with AB/2 electrode separation of up to 150m were used.
A resistivity distribution in the layers below 3m that satisfied (i.e.forward-calculated resistivities fell within one standard deviation of the corresponding field data set) the time lapse Wenner-Schlumberger data sets as well as the Schlumberger soundings from 1979 and 2007 was consequently used in all the starting models.

| RESULTS
The data analyzed in this work were collected between September 2012 and October 2015.An overview of the available raw timeseries, including gaps in acquisitions, is presented in Figure 5.

| Air temperature, precipitations, and ground thermal regime
The mean air temperature at the site in the analyzed period was is close to the active layer thickness on the same date, which was observed to be approximately 0.85m based on interpolation of temperature measurements from the MRC stick.

| In-situ unfrozen water content and freezethaw hysteresis
The unfrozen soil water content (Figure 5c) varies seasonally.The maximum volumetric water contents (up to 47% at 0.30m depth and up to 44% at 0.55m depth) are measured in the middle of June and again at the beginning of September every year.These two peaks are likely due to snow melt and a rainfall maximum, respectively-a general pattern common to all the three analyzed years.Between the two maxima and while the ground is above freezing point, increased water movement, evaporation, and precipitation events control the water content variation independently of the ground temperature.The minimum recorded liquid water content was approximately 10% at the upper, 0.3m sensor (corresponding to ground around temperature ca.À17.9 C), resp.approximately 7% at the lower, 0.55m sensor (corresponding to ground temperature ca.À14.3 C).
The zero-curtain effect sets in at the beginning of October and lasts for approximately 3-4 weeks at the depth of 0.30m.
The zero-curtain appears absent in the thawing branch of the unfrozen water content curve.This is because the length of the zerocurtain effect depends on the amount and availability of liquid water during the phase change process. 36,37nditions for lateral groundwater flow at the site exist only during the 4-month period between mid-June and end-October. 38In the timeseries of unfrozen water content (Figure 5c), at the end of each thawing event (mid-June at 0,3m depth), the soil briefly attains its highest unfrozen water content.This is likely due to melting of ice lenses that may initially limit the liquid water movement in the spring.
As more and more of the active layer thaws, fissures that occur in fine-grained soil following repeated cycles of freezing and thawing 30 cause higher hydraulic permeability and facilitate more rapid moisture redistribution in the thawing soil. 39,40In summer months June-August, drying and evapotranspiration lead to desaturation of the soil.This is confirmed by the observation that the water content decrease is smaller at the deeper soil moisture sensor.The end of August through September is the period of soil moisture recharge through rain, before the ground typically begins to freeze at the end of September.
The dependency of the unfrozen water content on temperature is distinctly different between periods of ground freezing and thawing.
Water content during freezing is consistently higher than during thawing at the same ground temperature (Figure 6a).The freeze-thaw water content hysteresis is obvious even during events of partial thawing in the frozen season (Figure 6b).In addition, freezing appears to be a slow process-likely due to higher quantity of available liquid water, freezing-induced hydrostatic and thus a more pronounced zero curtain effect-while thawing is comparatively fast.Nevertheless, the freezing and thawing patterns, respectively, remain the same every year (Figure 6b, c), and no noticeable inter-annual changes in total volumetric unfrozen water content were observed during the 3 years of our monitoring.
The freeze-thaw hysteresis is visually discernible in the temperature range from approximately À8 C up to the freezing point, even after decimating the 3-hourly unfrozen water content data to sampling frequencies of up to 5 days.However, to properly resolve the fast unfrozen water content changes during thawing, as well as for quantitative comparisons with the other soil parameters, at least daily logging frequency is needed.

| Seasonal enthalpy change at the site
To understand the significance of the observed unfrozen water content hysteresis for ground thermal calculations, we calculated the enthalpy change at the sensor location.Due to the observed hysteresis, enthalpy is not uniquely defined by the soil temperature but also depends on the unfrozen water and ice contents.We calculate the enthalpy change between two points in time, t 0 and t n according to the equation: where C eff ðtÞ is the effective heat capacity of the soil system at time t and given by: Here ΔH is the enthalpy change [J], T [K] is the ground temperature (from the soil moisture sensor, θ w [m 3 water /m 3 bulk ], θ i ½m 3 ice =m 3 bulk , and θ s ½m 3 grains =m 3 bulk are the volumetric contents of water, ice, and soil grains, respectively, when assuming fully saturated ground conditions.c w½J=m 3 water =K , c i ½J=m 3 ice =K, and c s ½J=m 3 grains =K are the specific volumetric heat capacities of water, ice, and soil grains, respectively.Table values were used (c w ¼ 4:190MJ=m 3 water =K, c i ¼ 1923MJ=m 3 ice =K, and c s ¼ 2:000MJ=m 3 grains =K) and assumed constant with temperature.L is the volumetric latent heat ð334MJ=m 3 Þ.T and θ w are measured quantities at given points in time (t), and θ i is calculated from the measured θ w using the equation: where η is the soil porosity [unitless] and T * [ C] is the effective freezing point of the soil (see Section 5.4).θ s is considered constant with time (1 À η).As a simplifying assumption, we neglect any volume change associated with phase change.
We chose the reference point of the calculations, t 0 , to represent the start of the hydrological year 2015 (2014-09-01), because it represents the approximate time where the active layer is fully developed, and thus the maximum enthalpy state of the system.The enthalpy change for any other point in time t n after t 0 can then be calculated by integrating along the path ðTðtÞ, θ w ðtÞ, θ i ðtÞÞ until time t n (in practice approximated by the discrete sum in Equation 2).We calculate ΔH for each day, at 18:00 UTC, of the full year following t 0 , and take the difference of the minimum and maximum values.We thus find the total yearly enthalpy change for the hydrological year 2014-2015 at the sensor location to be 174.7MJ.The observed unfrozen water content hysteresis translates to a hysteresis in the enthalpy change time series; the difference between the thawing and freezing branches was found to be 21.2MJ at À0.5 C.This corresponds to 12.1% of the observed yearly enthalpy change for the site.

| Parameterization of the unfrozen water content curves
Many algorithms and parameterizations of the unfrozen water content of frozen soils as function of temperature have been proposed (see 41 for a review).A simple and commonly used approximation takes the form of a power function θ w ¼ αjTj Àβ , [42][43][44][45] for T < T * < 0 C, where We assume the following relationship for the temperature dependence of the unfrozen water content: Here, θ w is the volumetric unfrozen water content of the bulk soil ½m 3 water =m 3 bulk , η is the porosity ½m 3 voids =m 3 bulk , ϕ is the volumetric unfrozen pore water fraction ½m 3 water =m 3 voids , S is the total water saturation ½m 3 water =m 3 voids in the unfrozen state, α and β are empirical constants describing the intrinsic freezing characteristics of the given soil, and T f is the freezing point of the pore water as a free substance.T * is the effective freezing point of the sample-the minimum temperature at which all the water in the sample is unfrozen ðϕ ¼ SÞ, as well as F I G U R E 6 (a) Field measurements (8 records/day, at 0.3 m depth) between 1 st September 2012 and 31 st August 2013 show the in-situ hysteretic variation of UWC during a course of 1 year.Desaturation due to water redistribution and drying is observed following complete ground thawing at the given depth (after ca. 25 th June).(b) The UWC follows the same freezing pattern year after year, as concluded from three consecutive freezing periods (1 st September to 28 th February, years 2012/2013 to 2014/2015).The freezing curve can be modeled with a set of soil-specific parameters ½a f , b f .Partial thawing is observed when the ground temperature fluctuates below 0 C. The slope of this partial thawing curve is the same as slope of the full thawing curve.(c) UWC during thawing (1 st March to 25 th June, years 2013-2015) during the three consecutive thawing periods can be modeled with a set of soil-specific thawing parameters ½a t , b t .
Equalizing these two equations, we obtain: We used Equations ( 5) and ( 6 The optimized parameters for the freezing and thawing seasons are summarized in Table 1; we observed that the α t and β t during thawing are ca.21% and 41% lower than α f and β f during freezing, respectively.Figures 6b and 6c show the volumetric unfrozen water content measurements from the three freeze-thaw seasons, as well as excellent agreement between the in-situ measurements and the modeled freezing and thawing curves.The two sets of α and β parameters calibrated on the freezing and thawing seasons in 2012/2013 predicted the volumetric unfrozen water content in the two years 2013/2014 and 2014/2015 within AE3%.

| Analysis of apparent resistivity
According to Hilbich et al., 7 analyzing apparent resistivity data before inverting them may improve interpretation of fast temporal changes of ground resistivity close to the surface.These fast changes otherwise tend to create inversion artifacts that complicate interpretation of results after inversion.The apparent resistivity variation in the winter season shows high sensitivity to sub-zero temperature and water content variations.
Measurements from the coldest period (February-mid-March) are often missing due to the extremely high grounding resistances restricting current transmission. 30During the thawing period, the increase in unfrozen water content is associated with a steady decrease in apparent resistivity, until temporarily plateauing in the month of July.In the following months of August and September the apparent resistivity further decreases until reaching its minimum in mid-September, thus corresponding well to the occurrence of the maximum active layer thickness as derived from borehole temperature data.

| Analysis of inverted resistivities
The timeseries of inverted resistivity soundings (Figure 7c) clearly show the processes of ground freezing and thawing.The freezing of the active layer during the months of October-November is reflected in a sharp resistivity increase, from ca. 100 Ωm to ca. 5 kΩm (Figure 8b).This is the period when the unfrozen water content decreases during the isothermal process of latent heat extraction, and the ground temperature records exhibit the zero-curtain effect (Figure 8a).During the coldest months of January-March, the highest resistivities reach the order of 10 4 Ωm in the topmost model layer.
The resistivity soundings successfully register several warmer events during the freezing period when the ground temperature is Unfrozen water content parameters determined by least-squares fitting of the θ w calculated using Equation ( 5) to the field-measured θ w from the Hydra Probe sensor at 0.3m depth.overall below 0 C.An example is the warmer period lasting for approximately 10 days in mid-November 2014 (Figure 7b, c).Resistivity soundings are thus very sensitive to small variations in unfrozen water content induced by sub-0 C temperature variations.They can therefore be successfully used to identify ongoing phase change in frozen soil close to its thawing temperature.
Between June and October each year, the inverted models suggest active layer resistivities in the order of 20-200 Ωm.The interpreted thickness of the active layer appears severely exaggerated though; the low resistivities associated with the unfrozen active layer are observed to depths larger than 2m in the inverted models, while according to the ground temperature data the thaw front does not exceed 0.9 m (Figure 7b).

| Sensitivity analysis of inverted resistivity model
To justify the use of the inverted resistivities in the comparison with ground temperatures and unfrozen water contents, we conducted a comprehensive analysis focusing mainly on the sensitivity of the inverted layer 2 resistivity (ρ 2 , depth range ≈ 0.15-0.32m) to variations in other model parameters.The analysis was conducted over one full year with high data availability (2014/15), and tested the effect of (i) changes to resistivities in the lower, fixed part of the model (layers 9-12), (ii) changes to the a priori and vertical constraints, (iii) fixed high resistivities in layers 6-9, and (iv) time lapse constraints.
For case 1, we found that ρ 2 was not sensitive to changes in resistivities of layers 9-12.These layers were tested sequentially with resistivities set to À50% and +100% of the values in the accepted model.The resulting mean relative differences in ρ 2 were less than 6%.Maximum relative differences up to 28% occurred during some periods of the year for changes to the layer 9 resistivity with little or no deterioration of the daily RMSE.The total RMSE (the RMSE calculated over all the stacked datapoints from all the daily acquisitions in the timeseries) change in all cases was less than 1%, relative to the accepted model.
Case 2 tested different a priori and vertical constraints.
Completely unconstrained inversion resulted in a slightly lower total RMSE (approx.À5%) but introduced unrealistic short-term variations in ρ 2 .These variations were suppressed with the introduction of vertical resistivity constraints in the upper eight layers of the model.Further tightening the vertical constraints (to STDF 1.5 and 1.2, down from 1.8 in the accepted model) significantly reduced the data fit in terms of both daily and total RMSE.
Case 3 consisted of four model runs.We imposed a maximum active layer thickness of 1m by fixing the resistivity of all the layers 6-8 to the same value: 3, 5, 10, and 20 kΩ m in the respective runs, while layer 9 was fixed at 1 kΩ m.This is not a physically accurate model as it does not allow any resistivity variation in layers 6-9, known to exist based on recorded yearly temperature amplitude, and the model fit is thus not optimal.The ρ 2 values were affected, with generally higher inverted resistivity values and larger day-to-day variation.However, the shape of the ρ-T relationship did not change significantly.
In case 4, we inverted the first, base acquisition using the initial model as per From our comprehensive sensitivity analysis, we conclude that ρ 2 -used in the quantitative temperature-resistivity-unfrozen water content comparisons-is strongly data constrained.Slight damping in the form of vertical constraints in the upper model layers is needed to avoid overfitting of the data.

| Relationships between inverted resistivity, ground temperature, and unfrozen water
To attempt a quantitative evaluation of relationships between inverted resistivity, ground temperature, and unfrozen water content, F I G U R E 8 Timeseries of inverted resistivity (solid black line) extracted from the second model layer between 0,15m and 0,32m depth ρ 2 compared with (a) the ground temperature (dashed gray line) and (b) the unfrozen water content (dashed gray line), both at 0,3m depth (averages of eight daily measurements).Zero-curtain periods are highlighted with the gray shading.The dependence of inverted resistivity on water content, (b), exhibits strong hysteresis effects, with resistivities during freezing one order of magnitude higher than during thawing at the same unfrozen water content (e.g.November 2014 vs. May 2015).
we compared the inverted layer 2 resistivity ρ 2 to the data from the upper moisture and temperature sensor at 0.3m depth.The petrophysical relationships at this depth are not influenced by the increased pore water salinity, which is only notable below 2m in the soil profile.
The rates of resistivity increase and unfrozen water content decrease were highest during the isothermal zero-curtain periods during freezing (Figure 8a,b).During these periods, the volumetric unfrozen water content decreased from ca. 44% down to ca. 28% in the active layer.Meanwhile, the resistivity increased from ca. 100Ω m to ca. 4kΩ m within the period of roughly 4 weeks.Resistivities at the lowest recorded ground temperatures (À17,9 C) and water contents (≈7%) could not be measured due to high electrode grounding resistances, with measured focus-one 46 grounding resistances reaching up to 170kΩ. 30However, the highest inverted resistivities were in the order of 10 4 Ω m at a ground temperature of ≈ À10 C. While there was a good correlation between fluctuations in inverted resistivity and unfrozen water content in the frozen period, the amplitudes of the resistivity fluctuations were proportionally larger than those of the unfrozen water content (Figure 8b).We speculate that variation in the electrode grounding resistances-correlated with surface temperature-contributed to larger amplitude of the inverted resistivity variations. 30ring thawing, the rate of decrease in ground resistivity was faster than the corresponding rate of increase in unfrozen water content (8c, f, i).In the temperature range between ca.À2 and 0 C and corresponding θ w range between 15% and ca.44% (i.e.throughout the major water content change), the resistivity observed during thawing was up to one order of magnitude lower than the resistivity during freezing (Figure 8b, e, h and 8c, f, i).These observations suggest a hysteresis in the resistivity vs. unfrozen water content and resistivity vs. temperature relationships.
In mid-June every year, when the frost table progresses below the soil moisture probe at 0.3m, the electrical resistivity further slightly decreases in spite of the decreasing water content.This may be an inversion artifact, as the depth to the high contrast layer boundary is increased.
Our comparisons revealed clearly defined relationships between ground temperature, soil moisture, and ground resistivity (Figure 9).unfrozen water content exhibits hysteresis too (Figure 9c, f, i), with inverted resistivities up to one order of magnitude higher during freezing than during thawing at the same unfrozen water content.

| DISCUSSION
Joint interpretation of the three types of comparatively dense monitoring data-unfrozen water content in the active layer, ground temperature, and electrical resistivity in the active layer and permafrost-offers new insight into permafrost and active layer processes.Naturally, more information also increases complexity and adds challenges to cross-interpretation of data collected at different temporal and spatial scales.We discuss some of the factors that could influence our measurement accuracies, potential effects on our interpretations, and the relationships between monitored parameters.

| Accuracy of unfrozen water content measurements
Dielectric permittivity is a well-established proxy for water content in soils.The calibration equation linking dielectric permittivity to water content (Equation 1) is based on the effective medium theory and on the optical path length of a single electromagnetic ray through the different medium components. 47It assumes that the geometrical shapes of the medium components are small relative to the wavelength, and evenly distributed at the scale of the measurement.The values of the two calibration constants depend on the relative permittivity of the individual medium components.Birchak et al. 47 developed the formulation for a two-component medium, but the formulation is extensible to multiple components.This implies three components (soil, water, and air) in an unfrozen soil and four components (with addition of ice) in a frozen soil.
The absolute accuracy of the in-situ measurements is difficult to quantify.Roth and Boike 48 estimated the accuracy of the θ w measurements by TDR to be within 0.02-0.05.Comparison with nuclear magnetic resonance (NMR) measurements in frozen soils showed that TDR 9,49 and FDR 9 devices are in good agreement with the NMR results.The dielectric moisture sensors are now commonly applied in in-situ active layer moisture monitoring protocols. 4,50

| Freeze-thaw hysteresis of unfrozen water content
Observations from both field (e.g. 4,8,502][53][54] ) agree on the presence of unfrozen water content hysteresis in soil freezethaw cycles.To examine the hysteresis in our field data, we first establish that the observed hysteresis is a real, measured phenomenon, rather than a result of influences by macroscopic effects or measurement errors.The factors that we consider are as follows: (1) effect of cryostructure development; (2) geometric effects and the sensor's volume of influence; and (3) heat conduction through the wires of the FDR sensor.
Cryostructure and development of ice lenses within the volume of influence of the FDR sensor could possibly result in a hysteresis effect.During ice lens formation, the ice lens grows from below.
Freezing of water at the lower boundary causes cryosuction, which brings additional water from below to the freezing front.The unfrozen water content below the ice lens will therefore be at, or below, saturation.In the thawing situation, with excess ice present in the soil, the melting of the ice lens will cause a situation with higher water content (than during freezing), until the excess water drains away.This hypothesized hysteresis effect would result in a pattern opposite to what we observe (higher water content during thawing than freezing) and would only occur while the ice lens is active (freezing or thawing).
It would therefore be observable at temperatures close to the freezing point, but not at lower temperatures.Therefore the cryostructure development does not explain the hysteresis at lower temperatures, nor during partial thawing.
The sensor itself could possibly introduce a hysteresis response, based on the fact that the measured temperature is a point measurement, whereas the permittivity measurement has a certain volume of influence within which the measurement is sensitive to variations in bulk soil composition.Periods when the freezing or thawing front moves through the volume of influence could potentially lead to an apparent temperature hysteresis which is not based on physical processes in the soil.This effect, however, would be limited to a narrow range of temperatures close to the freezing point, and thus unlikely to explain the consistent hysteresis observed at lower temperatures and during periods of partial thawing in overall frozen conditions.
It is theoretically possible that the temperature sensor in the FDR probe could be sensitive to heat conducted from the surface through its cable.This would lead to warmer temperatures measured during thawing than during freezing.In the freeze-up period (zero curtain observed), the records showed short periods of drastic air temperature increase, but no correlation was observed with the temperature measured by the sensor.We therefore conclude that the probe was not affected by such heat conduction effects.
Our field data showed consistent and repeatable patterns of the unfrozen water content hysteresis loop throughout the three monitoring years.Based on these considerations, we conclude that the observed unfrozen water content hysteresis is a true process and unlikely to result from any of the discussed error sources.
A number of explanations have been put forward for the freeze-thaw hysteresis phenomenon, including contributions by (i) metastable nucleation and effect of supercooling during freezing (both processes are absent during thawing), (ii) effect of electrolytes lowering the freezing point in the unfrozen water fraction, (iii) capillarity and its effect on lowering water chemical potential, (iv) adsorption effects (non-freezable bound water), (v) pore blocking (ink-bottle effect), and (vi) pore network interactions; contributions by a number of these processes are cited in multiple studies, for example, 51,[54][55][56][57][58][59][60][61][62] .Similar to the lab results of Tian et al., 54 we identified three zones of the hysteresis loop in our field data (Figure 6a); we discuss them together with the processes that are likely dominating the respective zones: • Zone 1-around the freezing/thawing point.Liquid moisture content above the freezing point is controlled by water movement in the soil and evapotranspiration rather than by soil temperature.
Right below the freezing point, the effects of metastable nucleation and supercooling dominate the hysteresis response. 54Zone 2-pronounced hysteresis.In the temperature range between the freezing (thawing) point down to ca.À8 C in our data, hysteresis is visually discernible with θ w during freezing consistently higher than during thawing at the same ground temperature.In laboratory studies, this θ w pattern was controlled by pore size and pore saturation: hysteresis was remarkable in larger pores (with smaller surface-to-volume ratio), while smaller pores (mostly containing bound water) showed practically no hysteresis. 58,61In terms of saturation, hysteresis was greatest in filled pores, while at low water contents-where most water is bound to pore surfaces-hysteresis effect was negligible. 59,61In natural, interconnected pore systems, pore network interactions contribute to the hysteresis through pore-blocking effects, as well as larger supercooling and consequently delayed ice formation in isolated, narrow pores. 56,58Zone 3-frozen soil.Below ca.À8 C at our site, the difference between freezing and thawing curves was not visually discernible, and θ w remained nearly constant with changing ground temperature.This is because at such low temperatures, unfrozen water mostly exists in the form of bound water, which does not contribute to the hysteresis effect. 59,63

| Resistivity imaging of the freeze-thaw processes
Resolution and accuracy of the inverted resistivity models are discussed in this section, as these have an impact on the quantitative simultaneous interpretations of resistivity, temperature, and unfrozen water content.
Resolution of the resistivity imaging.Resolution in the upper part of the resistivity model can be a reasonable concern.The resolution is determined by the minimum electrode spacing of 0.5m.The spacing was chosen from practical considerations of (i) the system specifications in terms of the number of electrode takeouts, and (ii) the tradeoff in terms of the electrode size-electrodes large enough to ensure good grounding during the freezing season yet small enough not to violate the inversion assumption of infinitesimally small electrodes, see considerations regarding the electrode size in (3).The optimal acquisition protocol should include the shortest possible electrode separations to resolve the near-surface variation, but the longest possible separation to accurately account for variations in the deeper layers.Furthermore, the depth sensitivity of a particular quadrupole configuration changes over the course of a year, due to the large variation and changing contrasts of resistivity at the site.Consequently, the depth of investigation (depth of maximum sensitivity) is shallower the deeper the active layer, as a larger part of the current passes through the active layer.Moreover, the high resistivity of the ice-rich permafrost layer decreases model sensitivity below this layer.Nevertheless, we argue that the consistency in the results of our sensitivity analysis (Section 5.7) documents that the system has sufficient resolution to quantify the ρ 2 variation, which is used for the comparisons with the unfrozen water content and temperature.
Resolution of the active layer thickness.Our inverted resistivity models indicate that the low resistivity associated with seasonal thawing extends to greater depths than the thaw depths known from the borehole temperature records.5][66][67] Very tight inversion constraints and a model specifically tailored to our particular data set was needed to accurately reproduce thaw depth propagation consistent with the available borehole temperature records.Such constraints would essentially prevent the model from capturing any future changes involving deepening of the active layer, and would therefore be of little use in general application.Rather, a critical choice and testing of inversion constraints in the deeper part of the resistivity profile (with little annual temperature variation) was our preferred solution to bring out realistic resistivity variations in the active layer.
Case 3 of the sensitivity analysis (Section 5.7) investigated the potential impact of the active layer thickness resolution on the ρ 2 values.
We concluded that the near-surface resistivity (ρ 2 ) was strongly data constrained as forcing a high resistivity boundary at 1m depthcorresponding approximately to the maximum thickness of the active layer known from the temperature records-did not significantly affect the shape of the ρ 2 -temperature relationship.Detection of the active layer boundary with resistivity method was systematically investigated by Herring and Lewkowicz, 68 concluding that maximum gradient-rather than a fixed threshold value-is a better indication of the boundary between unfrozen and frozen ground.

Resolution of resistivities in permafrost.
There was some ambiguity in resolving resistivity values in the parts of the model below the active layer.In summer time, the high-resistivity ice-rich layer at the top of permafrost was suppressed due to the effect of the lowresistivity layers above (the active layer) and below (the saline permafrost). 69In winter time, the effective depth of investigation of the longest Wenner-Schlumberger arrays in our data set was not deep enough to resolve the resistivity values of the saline permafrost layer.
Moreover, as 68 observed, resistivity of the frozen layers decreases model sensitivity at depth and increases uncertainty in the interpreted base of the frozen zone even when this lies within the depth of investigation.This effect is the most pronounced for laterally continuous, highly resistive permafrost, 68 just as is the case at our site.In addition, we speculate that the ground ice features in the top of the permafrost (Figure 2) may have caused pseudo-anisotropy, 70 while the top-down freezing and thawing patterns created a layered earth structure with strong resistivity contrasts, likely causing equivalency issues.Such effects can lead to over-estimating model depths and thicknesses.
These ambiguities can be addressed by (i) introducing a priori inversion constraints based on knowledge of the site's geology or ground resistivities, (ii) variable electrode spacing on the profile, with short spacing at the center (improving the resolution of the upper layers) and longer spacing at the extremities of the profile (improving depth of investigation and parameter resolution at depth), and/or (iii) adding more datapoints with different configurations (increasing depth sensitivity).
Our experience highlights the importance of critically interpreting inverted resistivity models from permafrost settings.Inverted resistivity models must be tested for equivalency issues before adopting a final interpretation.In this particular case, our main result-the shape of the ρ 2 -temperature relationship-was not affected and can be trusted, although the absolute values of ρ 2 may be offset (Section 5.7).
Time-lapse constraints vs. sequential inversions.Time-lapse inversion approaches (see review by Hayley et al. 71 ) are generally regarded as state-of-the-art methods for processing resistivity monitoring data sets.The benefit of these approaches, typically developed on data from temperate settings with smooth resistivity variations, has not been definitively established for acquisitions from permafrost settings 72 where ground presents high-contrasts resistivity environment.
Furthermore, resistivity acquisitions on permafrost are affected by errors that are more extreme and time-variable than those typically encountered in more temperate settings.This is notably due to changing quality of electrode grounding, measurement errors, and background noise 73 that result in different resolution, and in a need for different regularization of the individual acquisitions. 74,75Previous resistivity studies on permafrost 6,7,67 have not made use of the timelapse constraints and do not indicate whether these would have potentially improved the inversions.
Applying two types of time-lapse constraints (constraints on parameters and cascading initial models) in the inversion of our resistivity data set resulted in no benefit over the sequential inversion procedure that was adopted for the final interpretations (Section 5.7, case 4).Rather, the time-lapse constraints introduced significant smoothing of the inverted ρ 2 values, when the resistivity transition from high (frozen ground) to low (thawed ground) was delayed and extended in time due to the smoothing.This smoothing effect is particularly problematic for resistivity imaging of the thawing season, where the timeseries of unfrozen water content and ground temperature show that changes happen very quickly.Successful application of time-lapse inversion constraints relies on good knowledge and quantification of data error estimates and their variation throughout the time-lapse acquisitions.Large variations in electrode grounding resistances were observed in our data sets (previously reported in Tomaškovičová et al. 30 ), both in the temporal and spatial domains.
Introduction of data error weighting based on the grounding resistance measurements could lead to time-lapse inversion methods more suited for permafrost settings.
A reference model for quantifying the relative resistivity changes is not obvious to identify in a setting with naturally induced, cyclic yearly changes. 7A multilayered resistivity model is needed to capture the resistivity variation with depth, while the depth to the frost table would need to be constrained to address some of the resolution issues mentioned earlier.In high-contrast permafrost settings, a smoothing time-lapse constraint will make resolution of the freeze front boundary difficult.This could be addressed by implementing an adaptive a priori inversion constraint informed by an orthogonal timeseries, e.g.ground temperatures.

| Relationship between resistivity, ground temperature, and unfrozen water content
Time-lapse resistivity acquisitions have been successfully used to evaluate inter-annual changes in ground ice content. 6,7,14For such inter-annual comparisons, even irregular or comparatively sparse resistivity data sets have proved useful.With this study, we show that when calibration of the petro-physical relationships is sought, the dynamic and complex nature of the temperature-dependent processes necessitates sub-weekly acquisitions.
Apparent resistivities-instead of inverted resistivities-have been used to qualitatively illustrate the dependency between ground resistivity and temperature. 67Apparent resistivity values have the advantage of not being influenced by inherent inversion assumptions and errors.However, they are not quantitatively comparable to point observations at a specific depth (such as temperature and water content), as they correspond to volume averages over the total current path length.Thus calibration of quantitative relationships necessitates the use of inverted resistivities.We therefore stress the need to seek the best inverted resistivity models, as well as to critically weight potential geometrical influences of the measurement layout and effect of inversion process on the inverted resistivities.
Comparisons of indirect measurements from the ground surface (resistivity) with point observations at a specific depth (ground temperature and moisture content) inherently face scale and resolution discordance.This is because the inverted resistivity relates to a depth interval (thickness of the model layer, in our case ρ 2 depth range ≈ 0.15-0.32m),yet we treat it as point observation when quantitatively comparing it to the unfrozen water content and temperature measurements.This discordance is expected to be most prominent during the periods when the freezing or thawing front travels through the resistivity layer used for this quantitative comparison.
Similar to field observations of Doetsch et al. 10 and synthetic study of Herring and Lewkowicz, 68 no single resistivity value was associated with the freezing point-the resistivity increased approximately one order of magnitude during the isothermal zero-curtain process.The immediate increase in resistivity around the freezing point suggests that ion concentration in the pore water at the depth of ≈ 0.15-0.30m is small, and that water freezes relatively rapidly. 10ter the initial strong increase in resistivity at the freezing point, the further increase approximately followed a linear relationship between logarithm of resistivity and ground temperature.
During thawing, in the range between ca.À2 C and the freezing point T * , the inverted resistivity plotted consistently ca.one order of magnitude lower than during freezing against the same ground temperature and unfrozen water content values (Figure 9c, i).A similar relationship between resistivity and ground temperature was previously reported by Farzamian et al. 76 While studies on freeze-thaw electrical hysteresis are comparatively sparse, we exploit the similarity between drying/wetting vs. freezing/thawing processes to discuss possible explanations.
It has been recognized in laboratory conditions that because of hysteresis effects, soils may show markedly different physical properties at the same water content, depending on whether this water content was reached by wetting or by drying.This behavior was demonstrated for heat capacity by 77 and, 78 for thermal conductivity by 79 and for electrical resistivity by, 80,81 and recently by. 82In all the studies, the recorded thermal and electrical values were higher during drying (corresponding to freezing and lowering of conductive fluid content) than during wetting (corresponding to thawing and increase of the conductive fluid content).In field conditions, Overduin et al. 4 observed pronounced hysteresis of the ground thermal conductivity and unfrozen water content during three cycles of freezing and thawing.They also reported the sharp transition from frozen to thawed thermal conductivities during thaw, as opposed to a 3-month transition period during ground freezing.This is similar to the fast rate of thawing vs. slower freezing observed in our data (Figure 6).A recent review and discussion of hysteresis effects in sub-zero acoustic laboratory measurements are given by Lyu et al. 83 Laboratory results of 81 in particular-on electrical hysteresis of partially saturated sandstones-show striking similarities to our field experience with electrical hysteresis during freezing and thawing of clayey soil.Analogously, we identify three stages of the relationship between electrical resistivity and water content (Figure 9c): 1. Region 1-Frozen soil: Volumetric unfrozen water content up to ca. 14%, ground temperature up to ca.À2.5 C. Logarithm of resistivity is linearly dependent on ground temperature or unfrozen water content, no hysteresis is observed.
2. Region 2-Hysteresis region: Volumetric unfrozen water content between ca.14% and 44%.Hysteresis is pronounced in the freeze-thaw cycle, with resistivity during freezing consistently higher than resistivity during thawing at the same unfrozen water content.
3. Region 3-Thawed soil: Following complete ground thawing and up until the beginning of the new freezing cycle (i.e. when the ground temperatures are above the freezing point), changes in electrical resistivity are not primarily controlled by ground temperature variation.
We point out that the shape of freeze-thaw resistivity hysteresis we observed is contrary to the result we would expect for bulk resistivity estimates based on petro-physical relationships such as Archie's law 15 or geometric mean, 84 where higher conductive fluid content results in lower bulk resistivity.The fact that our freezing branch was more resistive suggests influence of processes that are not accounted for in these petro-physical models, possibly: 1. an inherent hysteresis effect in the physical processes governing the electrical current flow.This could be (i) due to different soil freezing vs. thawing patterns affecting pore-scale fluid distribution (suggested for the wetting/drying scenario by 81,85 ), or (ii) due to macroscopic scale thawing patterns, when during the thawing, a very thin but conductive surface layer develops on top of the still frozen soil, potentially acting as a preferential current flow path along the soil surface directly between electrodes 67 ; 2. surface conduction at the air/water interface resulting from charge density and zeta potential at the air/water interface 81,86,87 in a real, partially saturated (as opposed to assumed saturated) soil; 3. a system response, as a result of the measurement layout.Influence of the inversion procedure is unlikely, as the same hysteresis pattern emerged for the apparent resistivity vs. θ w relationship.
Furthermore, the magnitude of the resistivity variation in the coldest portion of the year is proportionally larger than the corresponding unfrozen water content variation, likely due to the added effect of grounding resistances.Once again, the commonly used empirical resistivity mixing relationships 15,16,84,88 do not fully capture the complexity of these real-world processes.
The phenomenon of resistivity variation with unfrozen water content in soils undergoing phase change thus warrants further investigation.Under in-situ conditions, naturally occurring freezing and thawing processes may lead to variation in pore sizes and pore-scale distribution of water, ice, and air.This implies that there may not be a simple relationship between the measured resistivity and in-situ unfrozen water content, as the relationship between these two parameters depends on the history of freezing cycles.Such observations suggest the need for site-specific calibrations of the resistivity-unfrozen water content relationship, as well as the need for systematic laboratory studies to support interpretations of freezethaw resistivity hysteresis.
Nevertheless, resistivity data collected from ground surface are often used in permafrost degradation quantifications and forecasts.It is therefore important to be aware that inverted resistivities differing by as much as one order of magnitude may be recorded at the same level of unfrozen water content due to freeze-thaw hysteresis.Large errors can be introduced into estimations of unfrozen water content from electrical resistivity if hysteresis effects are not considered.
This study is, to the best of our knowledge, the first to document such effects consistently in long-term field observations of resistivity and unfrozen water content.

| CONCLUSIONS
We described the results of a 3 years monitoring effort from highlatitude permafrost, comprising time-lapse acquisitions of ground electrical resistivity, unfrozen water content, and temperature.Our findings span from technical insights into operating autonomous resistivity acquisition systems to description of petro-physical relationships in soils undergoing phase change, and can be summarized as follows: • Autonomous monitoring station operation: Establishment of the remote control solution proved crucial for obtaining consistent timeseries of resistivity data, with acquisition gaps only during the coldest periods with poor electrode grounding (February-March).
These acquisition gaps were further reduced with optimized electrode design.
• Freeze-thaw hysteresis of unfrozen water content: A hysteresis pattern was identified in the seasonal variation in the unfrozen water content (θ w ), with θ w during freezing consistently higher than during thawing at the same ground temperature.An empirical power function linking θ w to ground temperature can be modified to account for this hysteresis when using separate parameterizations for the freezing (α f , β f ) and for the thawing parts (α t , β t ) of the θ w curve.
• Enthalpy hysteresis: The freeze-thaw hysteresis constitutes 12.1% of the yearly enthalpy change in the active layer.This is a nonnegligible proportion of the yearly enthalpy change for our site.
Therefore, separate parameterization of the freezing and thawing branches of the temperature-θ w relationship should be included in modeling of the thermal regime at the site.
• Critical interpretation of resistivity inversions in permafrost: Strong resistivity contrasts and ground ice features in permafrost landscapes may lead to overestimating resistivity model layer depths and thicknesses.Inverted resistivity models must therefore be interpreted critically, and tested for equivalency issues before adopting the final interpretation.Adapted inversion strategies are needed to manage equivalency issues in interpretation of resistivity data from cryospheric settings.
• Limitations of time-lapse inversion approaches for permafrost: Timelapse inversion strategies did not improve our inversion results compared to sequential (independent) inversions.Rather, timelapse constraints introduced unwanted temporal smoothing of inverted resistivities, reducing the resolution of the phase change processes of interest.To be effective in high-contrast permafrost settings with changing quality of electrode grounding, time-lapse inversions need to incorporate quantification of data error estimates, for example, based on grounding resistance measurements throughout the time-lapse acquisitions.
• Freeze-thaw hysteresis of ground resistivity: We observed consistently higher resistivity values during freezing than during thawing at the same ground temperature and unfrozen water content.While this observation was consistent with some other laboratory and field studies, it was contrary to the outcomes of empirical models relating bulk resistivity and conductive fluid content.The effect of the used inversion settings on the shape of the resistivity-temperature relationship was minimal and did not affect the strong hysteresis effect.Detailed laboratory measurements may establish more accurate quantitative relationships.However, we find the visualization of the resistivity vs. θ w relationship useful, because it is the field data that are typically used for modeling and predictions.This phenomenon warrants further investigations and awareness that resistivity differing by as much as one order of magnitude can be observed at the same θ w depending on the stage of the freezethaw cycle.
but ranges up to 2.3m in the general Ilulissat area (2016, personal communication from geotechnical site investigation for construction of the new Ilulissat airport).The upper 2-3m of the permafrost is very ice-rich, with thick and irregularly oriented ice lenses F I G U R E 1 Position of Ilulissat on the map of Greenland and location of the field site east of the Ilulissat airport.
such soils should have a freezing point of T f ≈ À 1:9 C. The measured ground temperatures at 4 m depth are in the range À3.0 to À3.4 C (see Figure 2c, and results in Section 5).

F I G U R E 2
Borehole information from various boreholes at the monitoring site: (a) Generalized soil column and classification.(b) Gravimetric water content (w) and porewater chloride concentration.The different water contents observed in boreholes ILU78020 and ILU2007-01 are primarily due to different sampling strategies.For ILU78020, the sample material was deliberately chosen to avoid excess ice, resulting in lower total gravimetric water contents.(c) Thermal regime visualized by minimum, mean, and maximum temperature observed by each sensor (trumpet plot) for hydrological years 2007/2008 (borehole ILU2007-01) and 2013/2014 (borehole ILU2013-01).The annotations are mean annual ground temperatures from the two time series at 4 and 6 m depth.(d) Example of ice content and structure in the transient layer near the top of the permafrost (ILU2013-01, depth 0.89-1.13m b.g.s.).
Overview of the components of the monitoring station in Ilulissat.(a) Winter view of the station toward northeast showing typical snow conditions.Ilulissat airport is to the left; position of the respective station components is indicated.(b) Summer view of the station toward north showing the vegetation cover.Airport facilities are visible in the background.(c) Inside view of the aluminum box housing the resistivity instrumentation, HOBO datalogger and Campbell datalogger (CR800 was upgraded to CR1000 in February 2014).White cables leading to electrodes of the resistivity array are connected inside the box to the electrode selector via two multi-pin connectors.(d) Cable system leading from the aluminum box to each electrode of the array, spaced 0.5 m.(e) Details of the soil moisture probe (Stevens HydraProbe) and the specific heat sensor (East 30).Image also shows the saturation conditions of the soil at the site-the digged up hole for installation of the sensors quickly fills up with water flowing in from the active layer.processingand interpretation of pseudo-1D data acquired using a Wenner-Schlumberger protocol in the time period between September 2012 and October 2015.The protocol was designed with the aim of achieving high resolution in the central part of the profile and of providing a good basis for comparison with the 1D borehole temperature data.To create the pseudo-1D protocol, we selected all the Wenner-Schlumberger configurations with a horizontal focus distance (horizontal location of maximum sensitivity, assuming a homogeneous half-space; for symmetrical configurations, this equals to the configuration center) in the central 4m of the profile (between 14.25m and 18.25m).This protocol design resulted in a total of 233 electrode configurations focused at 33 unique vertical focus depths (depth of maximum sensitivity of a given electrode configuration, assuming a homogeneous half-space) ranging from 0.26m to 5.69m and with four to nine individual measurements per each vertical focus depth.The injected currents ranged from 10mA in the thawed season down to 1mA in the frozen season.The current ranging was done automatically by the SAS1000.Each datapoint resulted from averaging two stacked measurements.Large amounts of resistivity data rely on automated filtering tools to detect and remove measurement errors.In our records, errors usually appeared as negative resistivity readings, high-value outliers, and measurements skipped due to high grounding resistances.We applied the first two steps of the filtering procedure described in Rosset et al.31 : (i) the technical filter-removing negative resistivity data points and measurements skipped due to high grounding resistance, and (ii) the magnitude filter-removing measurements higher than seven times the standard deviation of the technical-filtered data set.
1D modeling allows for comparability to the 1D borehole temperature data, and for control over the inversion process by constraining model geometry and resistivity variation in selected model layers using the prior geological information.Each day's 1D resistivity data set was inverted individually.The resistivity model domain consisted of 12 layers with layer boundaries representing the active layer, ice-rich permafrost, saline permafrost, and bedrock as inferred from the boreholes ILU78020 and ILU2007-01 (see Figure 2).The geometry of the resistivity model (the thickness of the model layers and depths to layer boundaries) was fixed throughout the inversion to reduce the number of inversion parameters.The upper 2 m were discretized into eight layers of logarithmically increasing thicknesses, to account for higher resistivity variation in the active layer and top of permafrost where most of the temperature changes happen.The remaining four layers corresponded to the lower part of the ice-rich permafrost, transition into saline permafrost, saline permafrost, and the bedrock (Figure 4).The initial model for each day's inversion had the same geometry and constraints but different initial resistivity distribution in the upper 2m.Resistivities in the model domain below 3m were fixed as little resistivity variation was expected (annual ground temperature amplitude < 1:2 C at 3m depth).This forced most of the resistivity variation to occur in the upper model layers-corresponding to the active layer and top of permafrostwhere also most of the temperature and phase changes occur.Initial resistivity values in the upper 2m resistivity model layers 1-8 were specified as the daily average apparent resistivities at the focus depths closest to the center of the given model layer.The resistivities in the upper 2m of the model were unconstrained with respect to the specified initial values, but resistivities between the neighboring model layers were loosely constrained in the vertical direction, to avoid overfitting the data (a standard deviation factor, STDF 33 , of 1.8 was used).

À2. 4 C
(min/max: À34.3 C / 19.1 C).Precipitations in the three observed years were between 235 and 265 mm (5a).August and September tend to be the rainiest months.First snow usually comes at the end of September or in October, and it typically melts throughout the month of May.The deepest ground temperature records are available from the borehole ILU2013-01, where at 6m depth the MAGT was À3.1 C (min/max: À3.1 C/À3.0 C between August 2013 and September 2014).Figure 2 shows a comparison between temperature envelopes in year 2007 and 2013.The MAGT at 4m depth increased by approximately 0.3 , from À3.4 (min/max: À3.8 C/À3.1 C, borehole ILU2007-01, between September 2007 and July 2008) to À3.1 (min/max: À3.2 C/À2.9 C, borehole ILU2013-01, between August 2013 and September 2014).Thaw depth probing along the ERT profile on the 10 th August 2013 showed thaw penetration down to 0.77 ± 0.06m (N ¼ 32).This F I G U R E 5 Timeseries of all the raw monitoring data (note the different depth scale for (c) and (d)).(a) Air temperatures, daily averages of 8 three-hourly records.Daily precipitations are presented with red bars.(b) Ground temperatures from the MRC probe, daily averages of 8 threehourly records at 16 depths.(c) Unfrozen water content from the two HydraProbe sensors placed at two depths in the active layer, daily averages of eight records.(d) Apparent resistivity from the Wenner-Schlumberger (pseudo-1D) configuration.Up to nine 4-electrode measurements at each focus pseudo-depth were averaged and are distinguished by the color gradient.The annotations indicate the major station servicing events-the rod electrodes replacement by meshes, and the addition of on-site computer with remote control, both contributing to better data completeness and quality.

α
and β are positively valued empirical constants, T [ C] is the soil temperature, and T * [ C] is the effective freezing point of the soil.
) to model the volumetric unfrozen water content and to provide a freezing curve parameterization reference useful for modeling.Due to the pronounced freeze-thaw hysteresis, two sets of parameters-½α f , β f for freezing and ½α t , β t for thawing-were used to realistically describe the water content variation.Based on the three available freeze-thaw cycles, we defined the freezing season as starting when the maximum active layer thickness is reached (beginning of September) and ending when generally the lowest ground temperatures are recorded throughout the soil profile (end of February).The thawing season spans the rest of the year, and only the days with ground temperature below the freezing point at a given depth are used for calibration.We used the freezing season September 1, 2012, to February 28, 2013, and the thawing season March 1, 2013, to August 31, 2013, to identify the corresponding sets of parameters ½α f , β f and ½α t , β t .The α and β parameters were calibrated by nonlinear least-squares fitting, with the cost function being the sum of squared differences between the measured and modeled unfrozen water contents.The ground temperature used in the θ w calculations comes from the temperature sensor of the FDR soil moisture probe.Although the soil water saturation varies at our site during the summer period, the soil was fully saturated at the onset of winter freezing and spring thaw each year in the study period.Thus for the purpose of the inversion, we consider S constant and unity in the entire winter season (September-February).In the summer season, only the data between March 1 and June 15 are used for calibration, as between mid-June and end August, the soil does not fulfill the saturation assumption.

Figure
Figure 5d shows timeseries of daily average apparent resistivities at each of the 33 vertical focus depths of the filtered Wenner-Schlumberger protocol.While the timeseries of just 3 years is too short to identify inter-annual tendencies or permafrost degradation, daily and seasonal changes clearly appear in the apparent resistivity data.The active layer thickness interval (the first meter of the soil profile) is represented by data points at seven vertical focus depths.The topmost focus depth is at 0.26m, which limits the resolution of resistivity changes right below the ground surface.

F
I G U R E 7 (a) Timeseries of air temperature, ground temperature (from the FDR soil moisture sensor) with unfrozen water content at 0.3 m depth.Zero-curtain periods are highlighted with the gray shading.(b) Timeseries of ground temperatures down to 1.5m depth, data from the MRC probe.(c) Timeseries of the inverted pseudo-1D resistivity soundings measured by the Wenner-Schlumberger protocol.The plot is zoomed to the upper 1.5m of the profile to highlight the resistivity changes in the active layer.(d) Timeseries of the apparent resistivity soundings.The apparent resistivity value for each day for every vertical focus pseudo-depth is averaged from up to nine acquisitions at that depth level.Date ticks indicate the first day of every month.

Figure 4 .
The inverted resistivities in the first eight layers of the base model were used as a priori values for the starting model in the following day, and constrained by a certain STDF during inversion.The approach was then cascaded into the remaining acquisitions in the data set, using always the previous day inverted model parameters as a priori values for the inversion.Loose time lapse constraints (STDF of 1.2) improved the total RMSE by less than 1%.Tighter time lapse constraints (STDF of 1.1, 1.05, and 1.01, respectively) caused deterioration of the total RMSE.All models with temporal constraints introduced significant smoothing of the ρ 2 values.The temporal constraints did not change the overall shape of the ρ-T relationship, neither did they improve the estimation of the active layer thickness.A time lapse inversion without temporal constraints while using the inversion results from a previous time step as a priori values for the following initial resistivity model gave nearly identical results to our accepted, independent sequential inversions runs.
The unfrozen water content hysteresis was observed in data from all three available freeze-thaw seasons 2012-2015 (Figure9a, d, g), and the inverted resistivity shows repeated hysteretic behavior in relation to temperature (Figure9b, e, h).This is expected, considering the temperature-water content hysteresis and resistivity dependence on water content.However, the relationship between resistivity and F I G U R E 9 Scatter plots depicting hysteretic nature of dependence of ground temperature, inverted resistivity, and unfrozen water content.Measurements of water content and temperature are from the HydraProbe at 0.3 m depth.The inverted resistivities are extracted from model layer between 0.15 and 0.32 m depth.Pictured are three complete freeze-thaw cycles (from 1.September to 31.August) 2014/2015 (a-c), 2013/2014 (d-f), and 2012/2013 (g-i); colors indicate the time scale.