Soil thermal properties during freeze–thaw dynamics as function of variable organic carbon and grain size distribution

Permafrost regions are experiencing increasing air temperatures, accelerating the thawing process, and thickening the active layer in summer. This can accelerate the release of greenhouse gasses into the atmosphere from the organic carbon stored in the permafrost. The long‐term thawing rates of permafrost below the active layer are governed by the soil thermal properties, the heat capacity, and thermal conductivity, which vary due to differences in grain sizes and distribution and organic matter content. Using nine column experiments comprised of fully saturated synthetic permafrost samples exposed to freeze–thaw cycles, the relative contributions of a range of soil grain sizes and organic matter contents on the soil thermal properties were investigated. The columns were subjected to a freeze and thaw cycle while soil temperatures were recorded in profiles. To infer the thermal properties from these experimental data, a numerical heat transfer model was used. The best fit between the observations and a batch of 5544 numerical models was used to find optimum values for permafrost thermal properties. The optimized heat capacity varied between 500 and 650 (J/m3 K) and thermal conductivity between 2.45 and 3.55 (W/m K). These optimized parameters were subsequently used to model a 100‐year permafrost active layer thaw scenario under warming air temperatures. Variations in the optimized thermal properties resulted in a time difference in thawing depth of 10–15 years and thawing depths varied between 9 and 10 m between the different optimized thermal properties at the end of the 100‐year scenario.


Vadose Zone Journal
carbon as greenhouse gasses CO 2 and CH 4 into the atmosphere and dissolution into groundwater (Mu et al., 2016;Schuur et al., 2008Schuur et al., , 2009Vonk et al., 2013). As a result of the thawing processes, the vast store of permafrost organic carbon is becoming part of the active global carbon cycle. Models used to investigate permafrost dynamics under different climatological conditions help to gain insight into thaw rates, thawing depths, and the reactivation of groundwater flows and associated carbon transport (Bense et al., 2009;Mohammed et al., 2021Mohammed et al., , 2022. However, parameterization of the thawing processes is challenging and requires judicious modeling choices (Lamontagne-Hallé et al., 2020). Few data are available on the actual key soil thermal properties, namely, the thermal conductivity and heat capacity. Parameterization schemes have been developed to calculate the thermal properties, but their performance is low due to limited data of soil thermal properties during the transition from frozen to thawed conditions, which is a subject that has received little attention (He et al., 2020)

Thermal properties
The thermal properties govern two heat transfer processes: heat transfer by conduction and advection. Conductive heat transfer is often the primary heat transport process in the permafrost subsurface, advective heat transfer is not activated as long as the groundwater is in a frozen state (Hayashi et al., 2007). Thermal properties consist of the thermal conductivity λ (W/m K) and heat capacity C (J/m 3 K). All constituents found in a permafrost soil, organic material, grain particles with a grain size distribution, and pore space partially or fully filled with soil moisture, have individual thermal properties, which can be averaged together in various ways to form the bulk or effective thermal conductivity λ v and heat capacity C v (Grenier et al., 2018;Midttømme & Roaldset, 1998;Yang et al., 2021). Effective thermal properties across the freeze-thaw interval are difficult to model, as the thermal properties are altered during phase change (Zhang et al., 2018). During this crucial phase, pore water becomes partially liquid, a process that is described using a soil freezing characteristic curve (SFCC) and varies for different soil types, saturation levels, and grain size (Kurylyk & Watanabe, 2013). The liquid porewater could activate groundwater movement and transport dissolved CO 2 and CH 4 . Heat transport processes during this phase are highly nonlinear and are strongly controlled by water saturation levels. There are still uncertainties on the link between soil properties such as grain size, organic matter content, and the resulting thermal properties during the phase change transition (He et al., 2020).
Research has been performed into the relation between the thermal properties and grain size, organic matter content, and

Core Ideas
• A numerical 1D heat transport model was used to simulate freeze-thaw transitions of permafrost column experiments. • The numerical heat transfer model accurately simulated laboratory freeze-thaw dynamics. • Numerically obtained thermal conductivity, heat capacity and porosity correspond well with field observations. • Variations in organic matter content affect both thaw depth and thaw timescales.
moisture content using the heating probe method in laboratory setups, resulting in models for the prediction of the thermal properties from such parameters (Abu-Hamdeh & Reeder, 2000;Dissanayaka et al., 2012;Hamdhan & Clarke, 2010;Zhang et al., 2018). Furthermore, Midttømme and Roaldset (1998) have investigated the effect of a range of various particle fractions on thermal conductivity and the thermal properties of mineral quartz fractions where they found a relation between grain size and thermal conductivity, and Mustamo et al. (2019) investigated thermal properties of peaty soils at the freeze-thaw interval.
Although much research has focused on permafrost thermal properties, these studies so far have not included the freeze-thaw transitions of fully saturated sandy soils. Here, we investigate permafrost built up out of sedimentary sand with a variability in grain size and organic matter contents over the freeze-thaw interval. To investigate these nonlinear heat transport processes, we used column experiments in a laboratory environment to collect data, and use a numerical heat transfer model to determine effective thermal properties.

Column experiments
Column experiments are typically used to recreate field conditions, where a soil column is radially insulated from its surroundings (Nagare et al., 2012b). The primary heat transport direction will be thereby vertical, creating a 1D heat transport system (Hayashi et al., 2007;Mohammed et al., 2014). Zhou et al. (2006) and Nagare et al. (2012b) stressed the importance of creating as best as possible 1D heat transfer conditions, and minimize radial ambient temperature interference. Previous column experiments by Nagare et al. (2012b) and Nagare et al. (2012a) collected a large active layer sample of Arctic peat and brought it to the lab to measure soil moisture and temperature profiles within a 1D soil column as it thawed. Similarly, Mohammed et al. (2014) used synthetically made peat samples to simulate active layer thaw using a 1D soil column experiment. The experiment by Watanabe et al. (2011) used unsaturated sandy, loam, and silt loam soils during the freeze-thaw experiments. In our research, we will focus on silty and sandy soils with variable organic matter contents under fully saturated conditions.

Research aim
The main aim of our study is to quantify the combined effects of various organic matter contents and soil particle size distributions on the bulk thermal conductivity and heat capacity at the freeze-thaw transition around 0 • C. This was achieved by creating a comprehensive representative data set for freezethaw characteristics in permafrost soils that highlights the differences in temperature propagation between samples consisting of contrasting grain sizes and organic matter contents. It is hypothesized that the thermal conductivity of the solid fraction would decrease with a higher organic matter content, due to the relative lower thermal conductivity of organic matter compared to sand, thereby changing the propagation speed of the freezing and thawing front. The samples were housed within soil columns and placed in a climate chamber, where temperature observations were gathered during simulated freezing and thawing events. This data set was then used to study the influence of these soil physical parameters on the thawing behavior and thus thermal properties of sandy permafrost with varying organic matter contents. The current study is a continuation of de Bruin et al. (2021), where the thermal properties of the active layer were determined using numerical models based upon solely temperature observations recorded at various depths over a 6-year period at the Quinghai-Tibetan Plateau (QTP). In this study, we use the same approach to determine the thermal conductivity and heat capacity using a modified version of the numerical heat transfer model by de Bruin et al. (2021). However, the temperature observations are collected using the 1D soil column experiments, instead of in the field. We aim to reconstruct the same temperature dynamics in the soil column experiments as observed in the field, but at a faster evolving timescale.

Experimental setup and instrumentation
A total of nine freeze-thaw experiments with synthetic soil columns were conducted within a temperature controlled climate chamber (Figure 1). The experiments were performed in two batches of five and four simultaneous soil column experiments. We use high-density polyethylene (HDPE) tubes to F I G U R E 1 Photo of two out of nine columns during the freeze/thaw experiment. Each column of 1.2 m high and 0.315 m diameter is wrapped in 0.12-m-thick insulation material held in place with duct tape and placed on a plastic pallet to isolate it from the ground. Also visible are the temperature sensor cables that come out of the insulation. house the soils, with an inner diameter of 0.315 m, a wall thickness of 0.028 m and total height of 1.2 m. Figure 2a shows the HDPE tube, which is on the sides and on the top insulated with 0.12-m-thick shell rockwool insulation. The passive insulation reduces steep radial temperature gradients, resulting in a small radial heat flux, which is monitored using temperature sensors in between the column and the insulation and accounted for during the analysis as it is part of the numerical model. At the bottom of the soil column, a detachable aluminium plate with a width and length of 0.4 × 0.4 m with a thickness of 0.02 m acts as a base, fixed with 12 heavy duty bolts and a rubber o-ring, sealing off the bottom. Simultaneously, the aluminium plate serves as a conductor to transfer heat to and from the base of the soil column to the surrounding air. In the climate chamber, the air temperature could be controlled from a range of −4˚C to 20˚C.
All soil columns were equipped with calibrated temperature sensors (Campbell Scientific T107). The sensors had a measurement range from −35˚C to 50˚C and accuracy of 0.01˚C. The sensors were positioned at heights of 0.05, 0.2, 0.35, 0.5, F I G U R E 2 (a) Cross-section of experimental setup with an inner high-density polyethylenecolumn equipped with (1) T107 temperature sensors, insulated by 0.12-m-thick shell insulation. (b) Translation of 3D soil column to the radial 2D heat transfer model domain. (c) Model domain and applied boundary conditions. Thermal properties assigned to materials can be found in Table 1. 0.65, 0.8, and 0.95 m above the column base in the middle of the sample. Sensors wires were fed through the HDPE column wall using watertight wirelock nuts. Temperature values were logged at 300 s intervals using a Campbell Sci. CR1000X logger. Two temperature sensors were placed in the climate chamber to measure the climate chamber air temperature.

Soil physical properties
Two freezing and thawing experiments were run with varying soil grain sizes and organic matter contents. The sedimentary grain sizes ranged from silty (D 50 of 175 μm) to medium sand (D 50 of 730 μm). To obtain these ranges, we used well-sorted quartz sand, which was sieved in calibrated sieves to achieve the three grain size ranges as reported in Table 1. Based upon the carbon contents reported in the studies by Siegert et al. (2002), Zimov et al. (2006), Strauss et al. (2013), and Dutta et al. (2006), we selected three average carbon contents for sedimentary deposits. We have added 1%, 5%, and 10% of the total dry sand weight in organic matter to the mixture. Table 1 gives an overview of the used soil grain size range and organic matter contents. Organic potting soil was selected as representative organic matter, composed from Swedish and Baltic peat T A B L E 1 Soil physical properties used in thisstudy.

No.
Particle grain size range (μm) with added perlite. The organic matter composition of these peat soils is considered representative for northern latitude peat permafrost soils. The nonorganic perlite is a part of the organic potting soil and has a low bulk density and high water capacity holding properties and was present in concentrations varying from 19 to 26 wgt% in five of the nine of the samples, and will contribute to the bulk thermal properties of the soil material. All samples were mixed in an industrial mixer for 1 h to create homogeneous samples, thereby preventing any measurements being affected by layering of the material. The columns were filled up to 1 m with soil material using the slurry packing technique with manual compaction (Lewis & Sjöstrom, 2010), which should result in a uniform packing and fully saturated samples, mimicking water saturated sedimentary deposits. Then, a layer of 0.1 m of water was added on top of the soil material. This provided a buffer medium on top of the sample that created a stable temperature boundary during the experiment. Furthermore, it had known thermal properties, which was required during subsequent modeling. After the freeze-thaw experiment, to determine the porosity and organic matter content, all nine soil columns were sampled at different depths using 100 cc sample rings. The porosity measurements helped to check the uniform packing of the samples and establish the water content within the samples. Each soil column was sampled with four sample rings, where two samples were collected at 0.8 m height and two at 0.5 m height within the column. Porosity was determined by measuring the mass difference between a fully saturated and fully dried sample divided by the volume of the ring. Drying of the sample took place in an oven at 105˚C for a duration of at least 24 h, evaporating the porewater but preserving the organic matter (Wilke, 2005). Organic matter content was determined using the loss on ignition (LOI) technique (Heiri et al., 2001;Wilke, 2005). Organic matter of a total of five samples is corrected for the perlite concentration, which was determined using the LOI technique.

Experimental temperature conditions
The temperature in the climate chamber followed an predetermined trajectory, mimicking temperature transitions at the freeze-thaw interval. At the start of the experiment, the climate chamber and samples had a uniform temperature of 3˚C. The temperature of the climate chamber was instantly lowered to −3˚C, starting the freezing process of the samples from the bottom-up with a vertically upward propagating freezing front. The climate chamber temperature was held stable at −3˚C until the samples were fully frozen and had obtained a uniform temperature distribution. After a stabilization period, the temperature in the climate chamber was instantly increased to 3˚C, starting the thawing process. This was a single step increase, which created a sharp temperature gradient, causing a quick melting process. The steep temperature gradient highlights any effects that the soil thermal properties have on the propagation of the freezing and thawing front as the processes occur over a shorter timescale. The freezer unit responsible for temperature management in the climate chamber periodically thawed itself to remove ice buildup on the freezer unit. This caused some temperature fluctuations in the air surrounding the soil column during freezing.

Permafrost heat transfer model
We use the software FlexPDE to solve the energy transfer equation considering heat conduction and phase change (https://www.pdesolutions.com) which is successfully applied in previous research (Bense et al., 2009;de Bruin et al., 2021).
The energy transport equation (Equation 1) uses the bulk thermal heat capacity of the various materials C v (J/m 3 K), and bulk thermal conductivity λ v (W/m K) to calculate the energy balance over the model domain. All materialspecific and model parameters can be found in Table 2. The model relies both on physical parameters that are coupled to material-specific properties such as the thermal conductivity and heat capacity of solid grains and water, and equations that describe system properties and behavior based upon empirically determined relations. In our study, bulk thermal conductivity (λ v ) and bulk heat capacity (C v ) were calculated by a volumetric average approximation of which the equations are listed in Table 2 (McKenzie et al., 2007), yet there are other methods as well to determine the bulk thermal properties. Here the individual thermal properties of water, ice, and solids form fractions of the total bulk thermal properties. This is based upon the porosity, saturation level, liquid water, and ice-fraction within the pores.
The SFCC is the curve that describes the fraction of pore water that partly freezes or thaws at subzero temperatures as a function of temperature, also known as the mushy zone (Kurylyk & Watanabe, 2013;McKenzie et al., 2007). In this model, SFCC is represented by the water saturation (S w ), which is an exponential temperature-dependent function, which is controlled by the dimensionless parameter W that determines the shape of the curve (Table 2). A small part of the pore water remains liquid at subzero temperatures, because of an increase in concentration of dissolved solids that are excluded from ice formation, and adsorbtive and capillary forces that depress the freezing temperature of the unfrozen water (McKenzie et al., 2007;Romanovsky & Osterkamp, 2000). This is set with the residual saturation fraction (S wres ). Parameter values for the SFCC are manually fitted once against the observations, and subsequently used for all the numerical models. Individual fitting of the SFCC parameter values would exponentially increase the amount of numerical models.
The model domain of the soil column was modeled as a 2D cylindrical half space (Figure 2b), thereby capturing heat transport in both the vertical and radial direction. This allowed us to better represent the actual heat transfer gradients in the column than would be possible with a 1D model, but limiting the computational load involved in a full 3D model T A B L E 2 Numerical model parameters (McKenzie et al., 2007), (Grenier et al., 2018 Figure 2c shows the model domain and applied boundary conditions. The left boundary condition represents the center rotation axis of the cylindrical domain with a zero flux boundary. All other sides of the model domain were forced by the smoothed observed air temperature. Temperature observation was recorded every 300 s and were smoothed to 2 h average temperatures. This was to dampen strong temperature fluctuations at the model boundaries and to reduce the computational load of the model. The entire freeze-thaw cycle was not modeled, but only the thawing interval starting at the moment that the air temperature was increased to 3˚C. This was because for the modeling and determination of the thermal properties, we use the undisturbed observed temperature signal as there are small temperature fluctuations in the climate chamber that occur only during the freezing process as the freezer unit periodically thaws.

Heat transfer model analysis
By varying the solid grain heat capacity (C s ), solid grain thermal conductivity (λ s ), and porosity (ϕ), a parameter space was created of a batch of 5544 unique combinations. Each model was forced with the 2 h averaged temperature observation. Subsequently, the model results were compared against the observed temperatures from the nine experiments. For each of the nine experiments, a subset of the total parameter space was made, selecting only the models with a porosity identical to the observed porosity, thereby optimizing the thermal conductivity and heat capacity, but constraining the total parameter space of the optimization based upon the observed porosity. Using the Kling-Gupta Efficiency (KGE) error measure (Equation 2), the best fit between the model and observed data was determined. This resulted in nine optimized parameter combinations for all nine experiments. The KGE error measure comprises of a bias (μ) between model (T sim ) and observation (T obs ), a scale misfit (σ), and the correlation (r) between model and observation (Knoben et al., 2019) and was successfully applied in previous numerical model optimization studies (de Bruin et al., 2021).

Optimum parameter scenario
To investigate the sensitivity of the thaw-depth rate to model uncertainty, we performed an analysis of the long-term effect of the parameter space on temperature propagation. The nine optimized parameter combinations based upon all parameter values were used to model a hypothetical 1D permafrost active layer thaw scenario over a 100-year period. The top boundary consisted of daily averaged air temperature observations observed over a 6-year period at the QTP, which were repeated over a 100-year period to create seasonal fluctuations. A linear temperature increase of 0.05˚C/year was added, which is representative for IPCC arctic climate warming scenario SSP5-8.5 (de Bruin et al., 2021;Constable et al., 2022). Initial temperature conditions were created using a steady-state model, which were subsequently used in the transient model as the initial conditions. The transient 1D model had a domain depth of 100 m where the bottom boundary consists of a fixed flux boundary of 0.065 W/m 2 , which was in range for geothermal heat fluxes at the (QTP) (Wu et al., 2010). The optimized thermal properties were applied over the entire depth domain, and optimized porosity over the top 20 m of the domain, a fixed porosity of 0.25 (−) was applied over the remainder of the domain as the optimized porosity was not representative for depths larger than 20 m. The resulting temperature models were analyzed to illustrate the development of the active layer depth and overall thaw depth over time.

Temperature observations
The temperature observations collected from the nine experiments in Figure 3 shows the temperature development at different depths within the soil columns during the freezing and thawing cycles over the two batches of five and four simultaneous soil columns. The observations in the figure were smoothed from 5 min observations to 2 h average temperatures. All observations within the soil columns show a similar temperature pattern. During the first phase I Figure 3b, the air temperature drops from +3˚C to −3˚C. The temperature in the soil column decreases to 0˚C and enters a phase change state. A number of temperature observations  show during initial freezing supercooling of the porewater for a short period of time indicated by the red circle in Figure 3a, where after temperatures increase to 0˚C (Akyurt et al., 2002;Ren & K.Vanapalli, 2019). During the propagation of the freezing front, temperatures remain stable at around 0˚C. During the phase change period, the soil column is partly frozen, and the mixture remains at a stable temperature until the pore water is crystallized. After crystallization of the porewater within the entire column (phase II) the temperature decreases to the set temperature applied at the bottom boundary. During phase III, temperatures within the column stabilize and respond quickly to external temperature fluctuations. During phase IV, the temperature in the climate chamber is increased to +3˚C and a thaw propagates upward through the soil column. During the thawing process, the observed temperature gradients are initially steep, and gradually decreases toward the end temperature of +3˚C. The observed air temperature shows significantly higher temperature fluctuations compared to the observation within the soil column, which are more buffered. When the thawing phase starts, the observed temperatures develop in a similar manner as during freezing. As the air temperature is increased to +3˚C, temperature within the column responds "quickly" and stabilizes at 0˚C withiñ 40 days. Pore ice starts to melt until all pore ice is melted. Thereafter, the 0˚C front propagates vertically through the column, with an initial steep temperature gradient followed by a more gradual temperature increase. During the freezing phase (I), the observations at 0.05 and 0.20 m depth show more fluctuations compared to the thawing phase (III), as a result of the periodic defrosting of the freezer unit. In Figure 3c, the temperature observation at 0.95 m depth suffered a malfunctioning sensor after 105 days.

Thermal properties
A numerical heat transfer model was run in two batches of 5544 models each. Resulting model temperature outputs have been analyzed using the KGE error measure to determine the best fit with the 2 h smoothed observed temperature data from the column experiments. Figure 4 shows the best fit between the observations and the optimized KGE model at a height of 0.35 and 0.80 m over the modeled thaw cycle. Furthermore, it shows how the thawing front propagates varying Panels e, f, and h show that there are small misfits in terms of the timing of the thaw front that a certain temperature is reached between day 50 and day 60. Figure 5 shows the optimized values for the thermal conductivity, heat capacity, and observed porosity determined by the numerical model and the organic matter content and grain size used during the experiments. There appears to be a decreasing trend in thermal conductivity with an increasing organic matter content. However, the thermal conductivity values with a D 50 of 730 μm show an increase in thermal conductivity. There is no clear relation visible between the thermal conductivity and grain size. The optimized values for the heat capacity fluctuate over a small range of values at the lower end of the modeled parameter space, with one value peaking at a higher heat capacity of 650 (J/m 3 K). There is no clear increasing or decreasing trend visible for the heat capacity in relation to either organic matter content or grain size. The observed porosity clearly shows an increasing trend in porosity with an increasing organic matter content. The highest porosity values are observed in the experiments with the largest grain size with a D 50 of 730 μm. The large grain size of D 730 μm has the highest thermal conductivity and porosity, but does not appear to have an effect on the heat capacity.

Porosity and organic matter
Each soil column experiment has been sampled at four locations to determine the average porosity and organic matter content. Table 3 shows the averaged porosity and organic matter content for each experiment alongside the numerically derived optimized thermal properties. The measured porosity values are all in the same order of magnitude and range expected under experimental conditions. There appears to be a relation between the measured porosity and particle grain size as can be seen in Figure 5. The porosity increases parallel for the different grain sizes.

Long-term model sensitivity
The nine optimized parameter combinations that best fit the observed data were used to model a hypothetical 100-year permafrost active layer thaw scenario. The results of the 100-year scenario highlights the long-term effect that different parameters have on permafrost thaw predictions. Figure 6 shows the interpolated active layer depths for the nine parameter combinations. During the initial 40 years, the active layer depth appears stable for all parameter combinations. However, after 40 years, the active layer depth increases at different rates for the various parameter combinations. A time lag appears for thaw depth for the different parameter combinations of 10 years at model year 70, and the time lag increases the following years. There is a spread in overall thawing depth after 100 years, ranging between 9 and 10 m. The parameter combinations with the largest time lag have the same value for the heat capacity (500 [J/m 3 K]) and identical values for thermal conductivity (3.35 [W/m K]), but differ for the porosity, which ranges between 0.4 and 0.5 (−), indicating that the model is very sensitive to variations in porosity.

Temperature observations
The temperature observations collected during the nine experiments ( Figure 3) illustrate clearly different freeze and thaw T A B L E 3 Numerically determined optimized thermal property parameter values with the lowest Kling-Gupta efficiency (KGE) error set out against the measured organic matter content and porosity.

Organic matter content (wt %)
Observed porosity (ϕ) (−) rates. This indicates that under near identical temperature and other experimental conditions, freeze and thaw rates are only influenced by the soil material and its thermal properties. Because the data are collected during two separate batches, small fluctuations in the air temperature make the temperature observations of both batches not directly comparable. Despite this fact, the usage of numerical modeling in both the vertical and radial direction allows us to use the observed air temperature as boundary condition and model the temperature development using a range of thermal properties. A crucial aspect of this research is the timing of the thaw events of all nine experiments. The delay in timing is a primary result of the difference in thermal properties. Figure 4 shows the difference in timing of the various temperature observations, which only differs a few days between the nine experiments. This is, considering the timescale of thawing processes (~70 days), a negligible difference. However, this difference in timing took place under a step temperature increase from −3 to +3. Under field conditions the air temperature transition will be less abrupt, which will increase the delay in timing between the modeled temperatures at the various depths. This is clearly visible in the 100-year future scenario (Figure 6), where the delay in timing is not in the order of days but up to 10 years.

Optimized thermal properties
The numerically determined thermal properties do not all show clear trends in relation to the organic matter content or grain size ( Figure 5). In the introduction, our hypothesis is that the thermal conductivity of the solid fraction would decrease with a higher organic matter content. This is due to the relative lower thermal conductivity of organic matter compared to plain sand. The numerically determined results show that there is indeed a decreasing trend in thermal conductivity with an increasing organic matter content. However, the trend is not very clear as the optimized values are scattered over a large range. Similarly, we would expect that the heat capacity increases with an increase in organic matter content due to the organic materials ability to absorb more energy per weight unit compared to quartsite sand. However, the optimized model results show little variations in optimal values and no clear increasing or decreasing trends. Next to that, the optimized heat capacity values are near the lower edge of the investigated parameter space with realistic values. de Bruin et al. (2021) showed that the heat capacity is a very insensitive model parameter compared to the thermal conductivity and porosity. Therefore, it is more difficult to numerically determine the optimum heat capacity, since variations in heat capacity only have a very limited effect on the temperature propagation. The observed porosity values show that porosity increases with organic matter content. As a result, there is a higher water content within the sample, which in turn increases the bulk thermal conductivity. The optimized parameter combinations determined using the KGE do somewhat coincide with the expected trends, but not convincingly. We have determined multiple factors that contribute to this. First of all, the sensitivity of the parameters varies greatly. A previous study by de Bruin et al. (2021) found that temperature is very insensitive to heat capacity variations, and thus difficult to determine accurately. The thermal conductivity and porosity are more sensitive and greatly affect temperature propagation. However, due to the dependency of thermal conductivity to porosity, there are multiple parameter combinations over a wide range that yield good model fits. Second, the optimization process used in this study uses a full 3D parameter space, which is evaluated using the KGE. This is based upon the concept that if there are multiple optimal parameter combinations possible, all possible combinations are explored. This is especially the case for correlated parameters. The limitation however is that multiple parameter combinations have good performing KGE fits, but the KGE only evaluates the error between model and observa- tion for each time step, but the timing between model and observation is also crucial. Because we evaluate the average KGE, the average KGE could be high, but there might still be a misfit in terms of timing between the model and observation at certain depths, as, for example, can be seen in Figure 4i. The optimization could be refined by taking a second error measure, such as the Russell's error (Russell, 1997) that evaluates the best fit between a model and observation based upon the time phase shift between the two signals (de Bruin et al., 2021), thereby evaluating the temperature values between model and observation using the KGE, and thereafter evaluating the timing of the phase shift between the model and observation.

Long-term effects of parameter variation on active layer depth
The temperature observations gathered during the experiments show only relatively small differences in temperature propagation and timing of the freezing and thawing front. However, the subsequent implementation of the optimized parameter values derived from the experiments into the hypothetical 100-year active layer model reveal significant variations in thawing depth and timing over a long timescale. During the initial 40 model years, active layer depth remains similar for all parameter combinations, as is shown in Figure 6. After 40 years, the active layer depth starts to develop at different rates. The two-parameter combinations that are spaced apart furthest have heat capacity values that are very similar to each other: 500 and 550 (J/m 3 K). The thermal conductivity and porosity values however vary. The spread in thermal conductivity along with porosity, are the primary contributors to the large spread in overall active layer thawing depth. This spread is clearly visible at, for example, model year 80, where the thawing depth between the two most outer model predictions varies between 5 and 7.5 m, which is a difference of 2.5 m.
The optimized parameter space appears to be quite "flat" when trying to determine the best fit for the experiments. There are multiple combinations that vary only slightly in their KGE performance. Based upon the subsequent 100-year future scenario we determined that the heat capacity has a low impact and thermal conductivity in combination with the porosity has a high impact on permafrost development on the long term.
The 100-year scenario ( Figure 6) highlights the implications of various thermal properties on the thermal development of the active layer, where especially the timing of thawing varies by over 10 years. In the 100-year scenario, the variation in timing impacts the permafrost environment on a larger scale as well (data not shown). Due to the increasing thawing depth of the active layer, a perennial thaw zone or intrapermafrost talik forms underneath the active layer (Walvoord et al., 2019). This would allow over time for advective heat transport to start next to conductive heat transport due to the activation of groundwater flow, potentially accelerating the development of taliks and ponds (McKenzie & Voss, 2013;Rowland et al., 2011;Sjöberg et al., 2016). This in turn creates new pathways for dissolved organic carbon stored in the permafrost to be transported to the surface where it can be decomposed and emitted as CO 2 and CH 4 (Ma et al., 2019;Turetsky et al., 2020).

Comparison of laboratory experiment with field observations
One of the main concerns with laboratory column experiments is the translation of laboratory scale processes to field scale processes, and the direct comparison of thermal properties found in the lab versus in the field. An important element of this translation is the SFCC, which is the curve that separates the solid and liquid ice fraction in the soil as a function of temperature (Grenier et al., 2018;Kurylyk & Watanabe, 2013;McKenzie et al., 2007). However, these curves vary for different soils and especially for experimentally packed soils in the lab (Devoie et al., 2022). To overcome this issue, we used a fixed SFCC that was adjusted to the observed lab data and subsequently applied during numerical modeling. By taking an SFCC that fits to the lab data, the modeled temperature propagation and resulting thermal properties well represent the observations.
In the previous study by de Bruin et al. (2021), an SFCC was used that was calibrated to the observed temperatures at the QTP. Because in both cases, the SFCC are adjusted to the specific conditions, the modeled temperatures fit to the observed temperatures. The resulting thermal properties found in both studies are therefore not impacted by a discrepancy caused by a mismatch in SFCC, and can therefore be compared independently of the SFCC.
Regarding the sensitivity, de Bruin et al. (2021) found similar results for model sensitivity to thermal properties. The optimized parameter values found in their study were based upon temperature observations from one location at the QTP. The study found optimal parameter ranges for thermal conductivity (λ s ) of 0.95 to 1.1 (W/m K), optimized heat capacity (C s ) of 875 to 1125 (J/m 3 K) and optimized porosity (ϕ) of 0.25 to 0.45 (−). In this study, we used nine temperature data sets instead of one. The optimal values found in the study by de Bruin et al. (2021) differ with optimal values found in this study for fine-grained sand with low organic matter content. This study found a higher λ s ranging from of 2.45 to 3.55 (W/m K) and lower C s varying between 500 and 650 (J/m 3 K). The main reason for this difference is the high porosity in this study compared to the study by de Bruin et al. (2021), affecting the optimized thermal conductivity. Site information of the QTP indicates silty sand soil with organic rich topsoil (Luo et al., 2018), which is in correspondence with the fine-grained sand.
Permafrost is subjected to increasing air temperatures, where the thermal properties highly influence the freezing and thawing rates as can be seen in the 100-year scenario ( Figure 6). Our results indicate a relation between thermal conductivity and organic matter content, where the thermal conductivity decreases with increasing organic matter content. This could imply that permafrost areas with a high organic matter content are thermally less conductive, compared to soils with a low carbon content, which would result in lower thawing rates for soils with a high organic matter content. However, caution should be taken when experimental results are translated to possible implications in the field, as the experiments only focus on a very specific conditions. The organic matter content range investigated in this study (0 to 4 wt%) is too limited to be able to extend this relationship to soils with higher organic matter contents. Despite this, large parts of mostly Yedoma permafrost areas contain low organic matter contents (Strauss et al., 2012). Future numerical permafrost modeling could focus on long-term effects of permafrost thaw in these Yedoma areas using the range in thermal properties found in this study. Especially the formation of perennial thaw zones (as was detected in the the 100-year scenario) indicates the reactivation of groundwater flows and associates dissolved carbons (Walvoord & Striegl, 2021). The variations in thermal properties affect the timing of the emergence of these groundwater flows. The timing is especially of importance since groundwater transports dissolved organic carbon that later on can be released as greenhouse gases CO 2 and CH 4 into the atmosphere.

CONCLUSION
In this study, we used laboratory soil column experiments as a method to investigate freeze and thaw processes around the 0˚C temperature range. Various soil grain sizes and organic matter contents produced significant differences in temperature propagation, indicating variations in effective soil thermal properties. With the help of batch numerical modeling, the temperature evolution of the column experiments was modeled and fitting parameter values were quantified. The optimized parameters all fall within the range of realistic values. The thermal conductivity and porosity were coupled parameters when evaluating the bulk thermal properties, and cannot be evaluated separately. The heat capacity was on the other hand an insensitive parameter compared to the thermal conductivity, making it difficult to accurately determine the precise value. Long-term model results of active layer development showed a time lag of over 10 years between different models. Furthermore, the overall thawing depth varied over 2.5 m at certain points during modeling. Both the time lag and thawing depth were the result of variations in thermal conductivity and the porosity. This indicated that these thermal properties were essential for long-term predictions and the associated timing of thawing events, and that it is thus crucial to consider the porosity and thermal conductivity. To conclude on the main research question, we were able to quantify the thermal properties of column experiments with various organic matter contents and grain sizes. Thermal conductivity appears to decrease with an increase in organic matter content, however there is still a high degree of uncertainty as only small organic matter contents were evaluated. Variations in heat capacity were not directly relatable to either soil grain size or organic matter content. This research demonstrated that we can simulate permafrost freeze-thaw dynamics in a laboratory setting, but need to consider the limitations regarding a direct comparison between laboratory and field data. Numerical models simulating permafrost freeze-thaw dynamics help to bridge the gap between the laboratory and the field and allow for comparison of permafrost thermal properties. Relations between thermal properties and soil physical properties help to increase our understanding of permafrost dynamics and future geohydrological conditions in cold climate regions, which is especially needed in a warming climate, where the size and timing of new emerging carbon fluxes from permafrost areas provide useful insights as they impact the ecology and society.

A C K N O W L E D G M E N T S
The project is funded by the Dutch Research Council (NWO) as part of an effort to investigate permafrost thaw and solute transport processes under grant ALWOP.237. We would like to thank the editor and reviewers for their comments, which greatly improved the paper.

AU T H O R C O N T R I B U T I O N S
Jelte de Bruin: Conceptualization; data curation; formal analysis; methodology; visualization; writing-original draft; writing-review and editing. Victor Bense: Conceptualization; funding acquisition; methodology; supervision; writing-review and editing. Martine van der Ploeg: Conceptualization; methodology; supervision; writing-review and editing.

C O N F L I C T O F I N T E R E S T S T A T E M E N T
The authors declare no conflicts of interest.