Stochastic representation of spatial variability in thaw depth in permafrost boreal forests

A simple stochastic representation of the spatial variability in thaw depth is proposed. Thaw depth distribution measured in the two larch‐type forests in eastern Siberia, Spasskaya Pad and Elgeeii, showed different spatial, seasonal, and interannual variability, respectively. Year‐to‐year variation in active‐layer thickness was minor in Spasskaya Pad compared to Elgeeii. A gamma distribution adequately represented both sites' thaw depth spatial variability as the cumulative probability. Thus, we developed a simple model using the gamma distribution that illustrates the spatial variability in thaw depth at any thawing stage as a function of a given mean thaw depth. A hierarchy of models was introduced that sequentially considered the constant state, linearity, and nonlinearity in the dependence of the rate parameter of the gamma distribution on the mean thaw depth. Although the requirements of the model levels differed between Spasskaya Pad and Elgeeii, the proposed model successfully represented the spatial variability in thaw depth at both sites during different thaw seasons.

7][8][9] An essential objective of monitoring the spatial and temporal variability in ALT was the determination of spatial representativeness. 4e spatial variability in thaw depth is a complex function of soil conditions (texture, components, water/ice content), vegetation, and organic layers.Thus, the deterministic model requires spatial distribution data on environmental parameters that are rarely available. 1For this reason, Anisimov et al. 1 proposed near-surface permafrost parameters, including ALT, as randomly spatially distributed variables consisting of both deterministic and stochastic components and developed a stochastic model to represent the ALT mean values and variances, assuming a normally distributed ALT.They showed that the ALT spatial variability measured at several sites in Alaska followed a normal distribution function.The distributions were not highly skewed, indicating that a normal distribution assumption of ALT was sufficient.However, Anisimov et al 1 also noted that the Shapiro-Wilk test for normality rejected the null hypothesis of normality in some instances.Therefore, it is uncertain whether a normal distribution adequately represents spatial thaw depth variability.Some thaw depth measurements showed skewed distributions with a long tail on the deeper side, particularly during the early thaw season (e.g., 2,10 ).Wright et al 10 pointed out that the logarithm of the thaw depths followed a normal distribution (i.e., log-normal distribution).Still, few studies stochastically represented thaw depth spatial variability for the early thaw season.Furthermore, because of soil surface constraints, for the probability distribution for the thaw depth at the shallowest limit, the normal distribution symmetric about the mean might fail to represent the thaw depth spatial variability.
The goal of this study was to represent spatial variability in thaw depth in a stochastic manner.Our study included manual thaw depth measurements at two boreal forest sites in eastern Siberia over several years at different warm-season times.We represent the observed thaw depth variability using a gamma distribution and propose a simple model to represent the spatial variability in thaw depth at any thawing stage as a function of a given mean thaw depth using the gamma distribution.

| Study sites and design
We measured the spatial distribution of thaw depth in two larchdominated forests in the middle part of the Lena Basin of the Republic of Sakha, Russia (Figure 1).The first area was the Spasskaya Pad Scientific Forest Station (62 ∘ 15 0 17 00 N, 129 ∘ 37 0 07 00 E, 214 m a.s.l. the Aldan River, approximately 300 km southeast of Yakutsk. 11The mean annual air temperature and precipitation observed at a nearby weather station (Yakutsk Meteorological Observatory) from 1981 to 2010 were À8.7 C and 236 mm yr À1 , respectively. 12om 2004 to 2008, this region endured approximately 1.5 to 2 times more precipitation than usual. 13High soil water conditions during this period adversely affected larch tree growth in Spasskaya Pad, damaging and killing some trees, 13,14 while no remarkable disturbance was observed in the forest at Elgeeii. 15 Iijima et al 14 found that damaged and subsequently dead trees in Spasskaya Pad were concentrated within a limited area of a "permafrost valley" with a deeper and oversaturated active layer, even in a small 50 m Â 50 m plot.For these reasons, we measured and compared the characteristics of spatial variability in thaw depth in these two forest sites in this study.
Cajander larch (Larix cajanderi Mayr) was the most dominant species at both sites, followed by silver birch (Betula pendula Roth.) and willow (Salix sp.). 16Partially, Spasskaya Pad consists of Siberian alder (Alnus viridis subsp.fruticosa (Rupr.)Nyman) 16 and Elgeeii consists of young Scots pine (Pinus sylvestris L.). 17 Both sites had similar forest floors that were dominated by cowberries (Vaccinium vitis-idaea L.) mixed with several herbs, such as red baneberries (Actaea erythrocarpa Small), and round-leaved wintergreen (Pyrola rotundifolia L.).[18] According to the previous soil surveys, the soils of Spasskaya Pad were classified as Cambic Cryosol (Calcic, Turbic) in the WRB 19 system 20,21 or Typic Haploturbels in the USDA soil taxonomy. 22,23On the contrary, the soils of Elgeeii were classified as Cambic Cryosol (Eutric, Loamic) in the WRB system. 24The humus horizon thickness did not exceed 5 cm on the Spasskaya Pad and averaged 10-15 cm in Elgeeii. 11ots of 50 m Â 50 m were set up at these sites (Figure 1).We routinely conducted multipoint thaw depth measurements at 25 points (the points of the closed circles in Figure 1) at both sites from 2016 to 2019 (with some exceptions; see Table 1).To capture more detailed spatial variability in thaw depth, we conducted thaw depth measurements at an extra 20 points on Spasskaya Pad (the points of open circles in Figure 1) in June 2017.

| Thaw depth measurements
We used a handheld dynamic cone penetrometer (TW-035, Sakatadenki Co., Ltd., Tokyo, Japan; hereafter, penetrometer) to minimize uncertainties in thaw depth measurements.The penetrometer consisted of a tip cone with a 60 angle and a 2.5 cm base diameter, a guide rod, a drive rod with scale, a knocking head, and a 5 kg slide hammer.The slide hammer, free-falling 50 cm along the guide rod, strikes the knocking head, which drives the cone into the soil.The advantage of this method is that it does not depend on the physical strength or skill of the measurer, unlike conventional measurements using a metal rod.Iijima et al 25 confirmed the applicability of a penetrometer to measure thaw depths by comparing them with traditional methods (metal rods, frost tubes, and soil temperature profiles) at three different sites in eastern Siberia.
In this study, we used the number of impacts required for 10 cm penetration N 10 as an indicator for determining the thaw depth.
where N is the number of impacts and Δd p (cm) is the corresponding increase in penetration depth.The procedure for measuring thaw depth was as follows: 1.The initial depth achieved by the penetrometer's weight was recorded as the initial value.

3.
Step 2 was repeated until Δd p was less than a given threshold ε.
4. When Δd p < ε, the number of impacts N was gradually increased, and the corresponding Δd p values were recorded.
5. The depth when N 10 reached 50 was defined as the thaw depth, D T (cm).
We set ε ¼ 1 cm in this study because the scale on the rod of the penetrometer was in 1 cm increments.
After removing the penetrometer, we inserted a rod with thermocouples into the existing hole and measured the vertical distribution of the soil temperature to determine whether the deepest point reached the frozen soil.Owing to broken penetrometer parts, we were forced to cease the measurements halfway.

| Analysis of spatiotemporal variability in ALT
The most straightforward way of analyzing spatiotemporal variability in ALT is to directly compare the measured ALT at each grid node over several years.This method shows the absolute interannual variation range of the measured ALT values.However, if the spatial mean ALT varies significantly annually, this may affect the interannual ALT variation range at each grid node.
To examine spatial variability in ALT at individual grid nodes over several years' time series, Hinkel and Nelson 6 proposed the normalized index of variability I v as follows: where Z avg (cm) is the spatial mean ALT for a particular year, and Z i (cm) is the node-specific value.Hinkel and Nelson 6 also defined interannual node variability (INV, presented as %) as the range in I v over several years, that is, the difference between the maximum and minimum values of I v in each node over several years.In addition, the grid-mean INV represents the average degree of variability in ALT over the entire recording period. 26According to previous results (e.g., 6 ), Smith et al 26

| Stochastic analysis of thaw depth spatial variability
In this study, we attempted to represent the spatial variability in thaw depth with a gamma distribution (e.g., 27 ).The advantage of the gamma distribution is that it can be represented by only two parameters, the shape parameter k (dimensionless) and the rate parameter λ (cm À1 ), and the mean and variance of the distribution are given by k=λ and k=λ 2 , respectively.Because the skewness of the gamma distribution is 2= ffiffi ffi k p > 0, the gamma distribution is positively skewed and converges with a normal distribution when k is large.The fitting of the gamma distribution to the observed thaw depth data was conducted using the R package "fitdistrplus" version 1.1-6. 28,29reover, because the number of measurement points is limited, the fitting of the function inevitably involves sampling uncertainties.
To assess these uncertainties, we conducted a nonparametric bootstrap analysis with 1,000 iterations and obtained the 95% confidence intervals (CIs) of k and λ of the gamma distribution.The R package "fitdistrplus" was also used for bootstrap analysis in determining CIs.

| Spatiotemporal variability in thaw depth
D T measured at each grid location in September, regarded as the ALT, showed consistent spatial variation in both Spasskaya Pad (Figure 2A) and Elgeeii (Figure 2B), irrespective of year.For example, the ALT at location C11 in Elgeeii was always shallower than that at other points (Figure 2B).This point was located in a depression slightly lower than the others, with high soil moisture and occasional waterlogging.Consequently, the higher ice content at this point than others necessitated greater latent heat to thaw, resulting in shallower ALT. 30 These results indicated that the thaw depths at individual grid points were forced by temperature and various local factors, and the point-specific ALT responded consistently across years, as suggested by Hinkel and Nelson. 6The consistent spatial variability in ALT over several years was also confirmed by the normalized index of variability I v (Figure 3).
The thaw depth D T measured at each grid location in Spasskaya Pad (A) and Elgeeii (B).
The year-to-year fluctuation range of ALT at each point was much smaller for Spasskaya Pad (mean: 5.7 cm, maximum: 13.5 cm) than for Elgeeii (mean: 15.2 cm, maximum: 35.0 cm) (Figure 2).The INV of Spasskaya Pad was also smaller than that of Elgeeii (Figure 3C).
This result implies that the conditions beneath the active layer differed between Spasskaya Pad and Elgeeii.We speculated that the transient layer 31 underneath the active layer in Spasskaya Pad constrained the maximum thaw depth.Despite such differences in the interannual ALT variability between the two sites, grid-mean INV was 3.7% for Spasskaya Pad and 8.2% for Elgeeii, both of which fell into "low variability." 26 contrast, D T variability during the middle of the thaw period poorly corresponded to ALT variability.Note that we measured D T near the grid points and were not precise at the same point every time, which would cause inevitable variability in measurements.
Nevertheless, considering that such uncertainty also occurs for ALT, this result indicates the processes determining the spatial variability in D T during the middle of the thaw periods might be much more complicated than that for ALT.For example, spatially varying snow disappearance dates may result in various spatial variability patterns in D T during the early thawing stage, which can make it different from that in ALT.

| Thaw depth and surface thawing index
The mean values of D T showed a clear linear relationship with the square root of the surface thawing index I TS (K d) in both Spasskaya Pad and Elgeeii (Figure 4; the numerical data are listed in Table S2).

| Fitting of the gamma distribution
Sorting the D T data for each measurement in ascending order represented the cumulative probability distribution, showing a sigmoidal shape (Figure 5).The distribution pattern differed at each measurement, but the distribution generally ranged wider in Elgeeii than in Spasskaya Pad and became wider when D T deepened.
To check the applicability of the gamma distribution to the mea- Section S2 for the detailed procedure).The cumulative distribution function (CDF) of gamma distribution was in good agreement with the cumulative probability of the merged D T data (Figure 6).The fitting of the gamma distribution was much better than that of the normal and Weibull distributions and similar to other asymmetric distributions (lognormal, Gumbel, and inverse Gaussian; see Figure S1 and Table S1).Although the fitting of the gamma distribution was not the best of these various distributions, we adopted the gamma distribution in this study because of the advantages mentioned in Section 2.4.
When bootstrapping was applied to the merged data in June  1.
was not always the case for the data of other measurements (see Supplementary Section S5).

| Modeling of spatial variability in thaw depth
The shape parameter k and rate parameter λ showed dependencies on the mean thaw depth D T (cm) for both the Spasskaya Pad and Elgeeii (Figure 9).Compared with k, λ is less variable against D T .Using this characteristic, we developed a simple empirical model to represent the spatial variability in thaw depth D T at any thawing stage as a function of the given mean thaw depth D T using the gamma distribution.
Because the mean of the gamma distribution is given by k=λ, k is expressed as the product of λ and D T as follows: Therefore, when D T is given (e.g., by the Stefan equation, see Section 3.2), we only need to parameterize λ to represent the spatial variability in the thaw depth, D T .We developed the following hierarchy of three-level models considering the various dependencies of λ on D T (i.e., constant state, linearity, and nonlinearity).
Model 2 considers the linearity of λ against D T .λ generally increased with D T in both Spasskaya Pad and Elgeeii (Figure 9) and is represented by a linear function as follows: In Elgeeii, Model 1 acceptably represented the spatial variability in D T at different thawing stages (Figure 10B).Although the cumulative probability of D T measured in Elgeeii had some variability in their sigmoidal shape, the distribution shape and its evolution with depth were reasonably reproduced by both Models 1 and 2. The difference between Models 1 and 2 was subtle; thus, Model 1 was considered sufficient for this site.This result indicates that Model 1 can be used as the first approximation for spatial variation in thaw depth at most sites where only the end-season thaw depth (ALT) was obtained.
However, in Spasskaya Pad, both Models 1-1 and 1-2 were insufficient to represent the spatial variability in D T across the thawing period, and Models 2 and 3 were required (Figure 10A).Model 1-1 where k 2:5 and k 97:5 are the 2.5th and 97.5th percentiles of k, λ 2:5 and λ 97:5 are the 2.5th and 97.5th percentiles of λ, and k obs and λ obs are k and λ obtained from the observed data, respectively.
Although the widths of CIs for k and λ varied significantly depending on the timing (mean thaw depth D T ) and site (see Figure 9), the relationship between NUR and sample size n showed similar characteristics regardless of site, timing, or whether k or λ, which falls along a single common curve (Figure 11; the numerical data are listed in Table S3).An exponential function of 1=n, obtained empirically from the relationship between NUR and 1=n, represented this curve.
The NUR with a sample size of n ¼ 25 had the highest number of measurements and varied more than other sample sizes, ranging from 1.204 to 1.396, including all NUR k and NUR λ in the Spasskaya Pad and Elgeeii.However, the mean and standard deviation were 1:331 AE 0:051, showing that most of the data concentrated within a narrow range, around the mean.Figure 11 and Equation (12)   show that the NUR increased as the sample size n decreased.The difference in NUR between n ¼ 25 and 45 was relatively small, whereas the NUR was significantly larger when n ¼ 17 and 18 compared to others.In other words, uncertainty did not decrease significantly with increasing sample size n when n ≥ 25, whereas it sharply increased with decreasing n when n < 25.This result confirmed that the sample size n ¼ 25 for regular measurements in this study was adequate.
The prediction function obtained in this study (Equation 12) can be applicable to other sites for uncertainty and sample size assessment.Note that the results in this study were obtained using

| Relevance and applicability of this study
We attempted to apply the gamma distribution to represent the spatial variability in thaw depth and to model the thaw depth spatial variability at any thawing stage.Our results revealed that the rate parameter λ in Spasskaya Pad and Elgeeii was low and close to each other in the early thawing stage but quite different at the end of the thaw season.The end-of-thaw-season λ was large in Spasskaya Pad, where some larch trees died due to high soil water conditions, while low and remained relatively constant throughout the thaw season in Elgeeii, where no tree damage was observed.These results indicate that the difference in λ at the end of the thaw season may reflect the different underground conditions.Such outcomes cannot be obtained if the normal distribution is adopted as in conventional studies (e.g., 1 ).
We assumed the presence of a transient layer in the soil of Spasskaya Pad.We need more verification with nearby borehole data or drilling for soil core sampling.Further discussion will be possible if our model is applied to the sites with available grid thaw depth data and borehole or soil core data.
Which probability distribution function to select for each site is still an open question.This study provided the applicability of the gamma distribution to the thaw depth spatial variability, but the data sets were limited to two boreal forests in similar environments.For example, Wright et al., 10 pointed out that the logarithm of the thaw depths followed a normal distribution (i.e., log-normal distribution) measured on a peat-covered permafrost slope in northern Canada.
Further investigation is necessary to discuss the appropriate probability distribution functions for various permafrost regions, including the requirements of statistical testing and plot scale issues (see Supplementary Section S4).Nevertheless, the gamma distribution and the model proposed in this study are expected to be strong candidates in their applicability and information obtained from k and λ.
; hereafter Spasskaya Pad), situated in a 200-year-old cowberry larch forest (Laricetum vacciniosum), located on a Pleistocene terrace on the western bank of the middle sections of the Lena River, approximately 20 km north of Yakutsk city.The second area was the Elgeeii Scientific Forest Station (60 ∘ 00 0 57 00 N, 133 ∘ 49 0 25 00 E, 202 m a.s.l.; hereafter Elgeeii) in a highly productive 180-year-old cowberry larch forest located in the third terrace of the left bank of the middle reaches of F I G U R E 1 Study site locations (left panel), measurement grids (middle and right upper panels), and photographs of forest floor conditions (middle and right lower panels) in Spasskaya Pad and Elgeeii.In the grid map, closed circles represent the regular thaw depth measurement points (25 points for each site) and open circles in Spasskaya Pad represent the additional measurement points in June 2017 (20 points).
presented a quantitative description of the mean INV as follows: (i) low variability for sites with the mean INV values of 0-19%, (ii) moderate variability for sites with a mean INV of 20-29%, and (iii) high variability for sites with a mean INV of 30% or more.

TheI
TS was calculated by the accumulated degree days of the daily mean soil temperature above freezing measured immediately below a soil surface (0 cm depth in Spasskaya Pad and 1 cm depth in Elgeeii) near the center of each plot.A single regression line through origin represented these relations well at each site, irrespective of the year, just as Anisimov et al 1 demonstrated.These results suggest that the simplified Stefan equation 32 can reasonably represent the mean thaw depth.
sured variability in D T , we prepared the merged D T data containing 45 data points in June 2017 on Spasskaya Pad from 25 regular data (June 15-16) and 20 additional data (June 22-23) (see Supplementary Normalized index of variability I v of ALT in Spasskaya Pad (A) and Elgeeii (B), and their interannual node variability (INV) in both sites (C).
2017, the obtained 95% CIs for k and λ were 14:43 ≤ k ≤ 35:87 and 0:229 ≤ λ ≤ 0:582, respectively (Figures7A and 7B).As a result, CIs around the CDF of the estimated gamma distribution were constructed (Figure7C) with a depth uncertainty of approximately 10-20 cm.The cumulative probability of the merged D T was within this uncertainty.The measured cumulative probability of D T at other times in Spasskaya Pad and Elgeeii was also mainly within the range of uncertainty (Figure 8).The range of uncertainty in Elgeeii was wider than that in Spasskaya Pad, partly because of the wider D T spatial variability in Elgeeii.Here, though the gamma distribution fitted the merged data in June 2017 better than the normal distribution, this F I G U R E 4 Relations between the square root of the surface thawing index I TS and the thaw depth D T in Spasskaya Pad (left) and Elgeeii (right).The symbols denote the mean values of measured D T , and the error bars represent the maximum and minimum values.F I G U R E 5 Cumulative probability distribution of the thaw depth D T measured at Spasskaya Pad (A) and Elgeeii (B).The boxplots shown together represent the distribution characteristics of the individual measurements, with the box showing the median and the 25th and 75th percentiles, the whiskers showing the 10th and 90th percentiles, and the cross showing the average.The dates shown are representative of each measurement period shown in Table

Model 1
provides λ as a constant.Because λ was less sensitive to D T in Elgeeii, we represented Model 1 for Elgeeii by the mean value of all measured λ. λ ¼ 0:292 ðModel 1for ElgeeiiÞ: ð4Þ However, in the Spasskaya Pad, λ increased significantly with D T .Large λ values at deep D T may have been caused by the effects of the transient layer.Moreover, a large λ value at D T ¼ 87:3 cm (July 2016) may have been possibly caused by the fewer data points (17 points).Note that the λ values measured in May and June were similar and close to the values in Elgeeii's Model 1 (Equation 4).In addition, the λ of June 2017 was the most reliable because it was evaluated from 45 data points and other λ from 25 or fewer data points.Therefore, if it is not for the transient layer, we expected that λ in May and June would represent all ranges of D T .Considering these circumstances, we tested two values for Model 1 on the Spasskaya Pad.Model 1-1 is the mean value of λ measured in May and June, and Model 1-2 is the mean value of all measured λ. λ ¼ 0:371 ðModel 1-1 for Spasskaya PadÞ: ð5Þ λ ¼ 0:959 ðModel 1-2 for Spasskaya PadÞ:

Model 3
considers the nonlinearity of λ against D T .The λ value in the Spasskaya Pad was significantly larger in September, whereas it was smaller in May and June.Therefore, the linear function cannot represent λ properly for the entire D T range.Furthermore, Model 2 (Equation 7) did not represent the most reliable λ of June 2017 evaluated from 45 data points.To represent nonlinearly changing λ, including this June 2017 value, we tested a nonlinear function in Spasskaya Pad.λ ¼ 0:238 exp 1:106 Â 10 À2 Á D T ðModel 3 for Spasskaya PadÞ: ð9Þ When fitting Model 3 to the measured λ, the λ in July 2016 was excluded because it had fewer data (17 points) and was considered less reliable.We did not test Model 3 for Elgeeii because k and λ in Elgeeii were satisfactorily represented by Models 1 and 2.

| DISCUSSION 4 . 1 |
Cumulative probability distribution of thaw depths D T measured in June 2017 on Spasskaya Pad.These 45 D T data were merged from regular (June 15-16, 2017, 25 points) and additional measurements (June 22-23, 2017, 20 points).The dashed line represents the cumulative distribution function of the gamma distribution fitted to the merged data.represented the distribution of D T in May and June reasonably but deviated from the results in September.However, Model 1-2 represented the distribution of D T measured in September but differed from the results in May and June.If λ is constant, the gamma distribution predicts a gradual increase in the variation range in D T with increasing mean thaw depth D T because k is proportional to D T (see Equation3) and the variance of the gamma distribution is given by k=λ 2 .However, in Spasskaya Pad, the range of variation in ALT was similar to that in D T during the mid-thawing season (Figure5A), suggesting that the ice-rich transient layer underneath the active layer restricted the maximum thaw depth.This explains the discrepancy between Model 1-1 (or 1-2) and the measured cumulative probability of D T (Figure10A) and why the change in λ should be considered in Spasskaya Pad.4 Effect of sample size on statisticsHow many data points are needed to capture the representative spatial variability in the thaw depth D T is an essential question for field researchers.If we could further increase the sample size, the reliability of the thaw depth spatial variability analysis would be further improved, but the measurement effort would also increase accordingly.In reality, the minimum sample size required to capture the representative spatial variability in thaw depth would be of interest.To assess the minimum sample size, we focused on CIs for k and λ, which we obtained by bootstrapping in Section 3.3.We defined the normalized uncertainty range (NUR) of a parameter as width of 95% CI divided by the value obtained from the observed data.The NURs of k and λ are given as follows:

7
Results of bootstrapping analysis (n ¼ 1,000) for the merged data in June 2017.(A) and (B) show the histogram of shape parameter k (a) and rate parameter λ (b) of the gamma distribution, respectively, and (c) shows the cumulative distribution function of the gamma distribution.Continuous and dashed lines represent the median and 95% confidence interval, respectively.F I G U R E 8 Results of bootstrapping analysis (n ¼ 1,000) for all regular measurements in Spasskaya Pad and Elgeeii other than June 2017 in Spasskaya Pad.The cumulative distribution functions of the gamma distribution are shown.Continuous and dashed lines represent the median and 95% confidence interval, respectively.nonparametric bootstrapping.If parametric bootstrapping is adopted, the obtained results can differ from ours.

F I G U R E 9 F I G U R E 1 1
The models of shape parameter k and rate parameter λ in Spasskaya Pad and Elgeeii as functions of the mean thaw depth, D T .The measured values, bootstrapping medians, and 95% confidence intervals of k and λ are also shown.To simply represent the spatial variability in thaw depth in permafrost regions, this study discussed the applicability of a gamma distribution to the measured thaw depth distributions in two larch forests in eastern Siberia.The thaw depth spatial variability characteristics differed between Spasskaya Pad and Elgeeii, with less variation in Spasskaya Pad, particularly for its seasonal maximum (i.e., ALT).In Spasskaya Pad, a transient layer underneath the active layer is speculated to constrain the maximum thaw depth.The gamma distribution well represented the measured thaw depth spatial variability at both sites with a 95% CI.We found that the shape parameter k and rate parameter λ of the gamma distribution depended on the mean thaw depth.Based on this finding, we developed a simple stochastic model that uses the gamma distribution to represent the spatial variability in thaw depth at any thawing stage as a function of the given mean thaw depth.This model consists of three-level models expressing λ dependency on the mean thaw depth.Model 1 represents λ by a constant, Model 2 considers the linearity in λ, and Model 3 considers the nonlinearity in λ.Although the requirements of the model levels differed between the Spasskaya Pad and Elgeeii, the proposed model successfully represented the spatial variability in thaw depth in both sites at different thaw seasons.This case study was conducted in two boreal forests in eastern Siberia within 350 km of one another, measured on a limited plot scale of 50 m Â 50 m.Therefore, further investigation is required to discuss the applicability of the gamma distribution and model proposed in this study to other boreal forest regions and regions other than boreal forests, especially tundra, and to confirm the spatial variability in larger areas (see Supplementary Section S4).Moreover, our model's coefficients for the rate parameter λ are expected to be represented by other environmental conditions, such as climate zone, soil types, and plant functional types.This may cultivate a further understanding of phenomena and allow robust modeling of active-layer dynamics and their impact on various processes in permafrost regions under changing climate.F I G U R E 1 0 Examples of models of the cumulative distribution function (CDF) and probability density function (PDF) that represent the spatial variability in thaw depth D T at different timings during the thaw season.Sample size dependency of the normalized uncertainty range (NUR) for shape parameter k (NUR k ) and rate parameter λ (NUR λ ) in Spasskaya Pad and Elgeeii.The dashed line shows the common curve fitted to all the NUR data represented by an exponential function of the reciprocal of sample size n.