RESULTS AND DISCUSSION

Model initialization

The minimum number of spin-up cycles needed to reach stabilization for soil temperature and moisture content under different initial conditions is summarized in Fig. 5 . Several points can be observed:
Both sites showed a noticeably larger thermal memory compared to the hydraulic one, as found by Elshamy et al. (2020). It worth noting that the WH site did not form permafrost under warm conditions (i.e. soil layers temperatures > 0°C), as discussed later in this section.Possible Position of Fig. 5 .
The behaviour across soil layers varies as the spinning-up process continues (Fig. 6 ), which is attributed to the function of each group of layers in water and energy exchange, including the root zone extent, the percentage of organic matter, the depth of organic soil (denoted ODEP thereafter), SDEP, and most importantly, the memory of the system (Taylor et al. , 2006). In general, temperature profiles change rapidly within the first few cycles and then stabilize with no significant changes after 1000 cycles, or less, depending on the initial water content and the climate of the spin-up year. After stabilization, most simulations showed further minor fluctuations, which can be attributed to model numerical instabilities, but their impact on the simulations is minimal. Soil moisture profiles depict the same behaviour with a relatively higher variation for the few first cycles. For example, Fig. 6A shows the pattern of change of each soil layer’s temperature, water content, and ice content over the spin-up cycles at the JMC site, under dry spin-up year climate conditions, and relatively dry initial soil content (Exp 7: 25% liquid + 25% Ice). The fully organic soil layers (i.e . layer 1:8) exhibited a sharp change within the first 50 cycles, followed by relatively insignificant variations for soil temperature and water contents. The lower soil layers down to SDEP (i.e . layers 9:14) had considerable temperature variations. In contrast, only the layers at the interfaces at ODEP and SDEP displayed significant oscillations in soil water and frozen content, due to the sudden change in soil properties that led to numerical artefacts – SDEP was not nudged/relocated to the nearest layer’s boundary in the current study. Soil layers below SDEP (i.e . layers 15:25) showed a diminishing pattern of variations with depth, with no variations at the lower boundary where the geothermal flux is set to zero.
Fig. 6B presents the initialization results at the WH site under wet spin-up climate conditions and relatively dry soil conditions (Exp 17: 12.5% liquid + 12.5% Ice). The MSO layers (i.e . layer 1:6) stabilized after a few cycles (~10) for temperature and water content, with negligible variations afterwards. As for the JMC site, lower soil layers down to the SDEP (i.e . layers 7:16) showed the main fluctuations for temperature (stabilizes after 100 cycles), while the interfaces at SDEP and ODEP drive the variations for liquid and frozen water contents. Another example for the WH site is provided in Fig. 6C , which shows the spin-up convergence under average climatic conditions and high initial ice content (Exp 10: 0% liquid + 75% Ice), noting that this experiment did not form/initialize permafrost at the end of spin-up. This is mainly attributed to the characteristics of the spinning year’s climate condition (i.e. average), rather than the specified initial soil moisture – all simulations under average and wet climate conditions did not form permafrost (discussed later in this section).
Fig. 7 summarizes the required spin-up cycles needed by the three state variables to reach equilibrium simultaneously for each climate year and each initial soil storage scenario. The two sites required 200–1000 spin-up cycles to perform an appropriate model initialization, depending on the spin-up year’s climate, as shown inFig. 7A . However, employing average and warm climatic conditions at WH site led to faster stabilization (< 200 cycles) without forming permafrost. Initializing the two sites under wet climate resulted in the shortest spin-up needed to form permafrost, compared to the other climatic conditions. On the other hand, grouping the initialization results regarding initial soil storage provides further insight into the impact of saturation level and the associated partitioning into liquid and frozen contents ( Fig. 7B ). Extreme saturation conditions, i.e. 100%, 75% and 0%, yielded the longest required spin-up cycles among all configurations at both sites. Average initial total soil moisture, i.e. 50% and 25%, required relatively shorter spin-up to stabilize, with best performance for 25% total soil water content that corresponds to field capacity condition. Further, Exp 16 (25% liquid + 0% Ice) and Exp 18 (18.75% liquid + 6.25% Ice) showed the minimal spin-up needed to initialize and form permafrost at the two sites under different climatic conditions. Note that the partitioning of liquid and frozen contents played a central role in the spin-up rather than the overall degree of saturation, as depicted in the first five experiments (100% saturation with different partitioning).
As mentioned earlier, the annual average air temperature and annual total precipitation were utilized to distinguish the various climate conditions (Sapriza-Azuri et al. , 2018). However, employing such mean/total measures could be insufficient and misleading, ending up in an unrealistic simulation. Fig. 8 shows the daily-evolution of air temperature and ‘cumulative’ precipitation at the WH site. Our selected ‘average’ year (i.e . based on annual total/mean) exhibits wet winter precipitation, while the fall’s air temperature included some relatively warm events. Further, the ‘average’ year had a cold-dry winter and spring, which led to relatively lower accumulation of snow on the ground. Likewise, the representative warm year recorded the lowest air temperature at the onset of winter and spring, coincident with relatively the wet- and dry-year conditions, respectively. Since WH site was configured without heavy organic soil (i.e. peat) and the surface vegetation was shrub canopy, the presence of sufficient snow is the major element for buffering/regulating the impact of external forcing on soil layers, and hence, the formation of permafrost (Dobinski, 2011). Another example can be seen in the representative ‘wet-year’ having the lowest precipitation volume for half of the year, coincident with the selected warm-year condition at the same period. Such intra-annual variations could affect the initializations drastically, as for the case of the wet year that was a combination of dry and warm conditions in winter and fall, disallowing any insulation by snow and providing extra heat flux at the vegetation-soil interface.
It worth noting that for the case of WH, our selected spin-up year with ‘average’ climate is the second coldest year among the five years (seeTable 4 ), with an annual average air temperature of -2.2°C did not form permafrost, while the other warmer years, i.e.designated as Wet (-1.14°C) and Dry (-1.19°C), have successfully formed permafrost after a longer spinning effort. Fig. 9 compares the daily evolution of the external forcings (i.e. precipitation, air temperature and accumulated snow depth) and the associated response of the soil state-variables (i.e. soil temperature and moisture contents) under average climate (left panel) and cold climate (right panel) for the same initial moisture condition (Exp 1: 100% liquid + 0% Ice) at the WH site. The peculiarity of the last week (or ten days) of September is the main reason for not forming permafrost under the ‘average’ year. In detail, the air temperature rose (~10°C) and was accompanied by a rainfall event (~20mm), which seeped down to the soil system and warmed it up – liquid content increased and hence the system’s thermal conductivity. Subsequently, the warming propagates for few spin-up cycles until the system reaches a self-consistent state. The spin-up needed to attain stabilization (without forming permanent ice) is controlled by the specified initial moisture conditions, which vary between 10 and 200 cycles for both the average and warm initial climates for WH setup.
Possible Position of Fig. 6 .
Possible Position of Fig. 7 .
Possible Position of Fig. 8 .
Possible Position of Fig. 9 .

Uncertainty propagation

In this section, we propagate the uncertainty of initialized model states (at the end of spin-up) through the simulation period and quantify the resulting uncertainty using two error metrics describing the simulation quality. We also focus on the temporal variation of different simulated aspects of permafrost. Fig. 10 shows the envelopes of soil temperature, liquid content, and frozen content profiles at the end of spin-up for the two sites, from all initialization experiments. Given that the presented results correspond to the last day of the spin-up “Sep 30th”, the following points can be noted:
Fig. 11 shows the impact of spin-up conditions on two performance metrics of the simulated permafrost, noting that the length of available record varies between the two sites. The bias in ALT (ɛALT) is insensitive to the climate condition of the spin-up year for the JMC site –ALT is overestimated by ~ 0.3 m, except for the cooler condition, which underestimates ALT by ~ 0.3 m (see Fig. 11 ). Initialization with a dry soil (i.e. Exp 21: 0% Water + 0% Ice) caused an apparent overestimation of ALT of the order of 0.8 – 1.2m. On the other hand, the simulated ALT at the WH site exhibits a general tendency to underestimation, by 0.2 - 0.8m, underlining partially the insufficient insulation provided by the MSO. It is worth noting that the experiments that failed to establish permafrost during the spin-up procedure developed it within the simulation period. The average- and wet-based experiments overestimate ALT by 0.25m, unlike the rest of the experiments that underestimate it by a similar magnitude (see Fig. 11 ). Experiments with the fully dry soil condition (i.e . Exp 21) did not show any differences between sites, which can be associated with the role of soil-texture and land-cover parameterizations.
The second performance metric used is the mean RMSE of the annual temperature envelopes (see Fig. 11 ) - calculated for the entire profile over every year, whenever observations are available. Since the depths of observations are not the same as used to configure the soil column, the RMSE is calculated by interpolating the simulated temperature envelope to the observation depths using the nearest neighbour method. While for JMC, ALT is overestimated with higher RMSE of the annual temperature envelope, varying by ~ 2.0°C at each simulation year, the WH setup, which tends to underestimate ALT, provides more accurate annual envelopes with a range of variability around 0.5°C. Fig. 12 shows the annual RMSE of Tmax and Tmin at the two sites. Again, the two sites tend to have cooler Tmin envelopes than observed (not shown ), resulting in higher errors for the cold part of the mean temperature profile. Regarding the annual maximum temperatures at JMC, most simulations produce the warm envelope with lower RMSE, with an upper limit of the error of +1.5°C. The annual minimum envelope at the WH site shows the opposite behaviour as the mean of simulations is closer to the upper limit of variability for all the experiments, with a narrower magnitude compared to the JMC site. These two examples highlight how the uncertainty inherent in the initialization of LSM can change the quality of the simulated soil temperature envelope.
Possible Position of Fig. 11 .
Possible Position of Fig. 12 .
Fig. 13A summarizes the simulated MAGTp calculated at the top of permafrost (bottom of the active layer). JMC showed cooler MAGTp compared to WH. The impact of the initial soil moisture is not identical among different climate conditions at the two sites, underscoring the influence of the driving climate on MAGTp profiles and the propagation of its uncertainty. Cold conditions at JMC tend to produce cooler permafrost (cooler MAGTp); a similar observation can be noticed for cold and dry climates at the WH site. The annual range of uncertainty of MAGTp is ~ -2.5°C at both sites.
Analyzing the offset effects (i.e. thermal and surface) reveals the high impact of the thermal offset occurring within the active layer. Initial soil moisture and climate conditions affect the cooling/heating flux to the permafrost, ranging between 2-3°C annually among all the experiments at the two sites. On the other hand, the ‘surface offset’ above ground due to the accumulated snow in winter and the standing canopy in summer shows insensitivity to the initialization condition (i.e. initial moisture content) with an insignificant range of variability among all tests (figures not shown ).
The day of maximum thaw (ALT-DOY) is among the most sensitive variables to the model’s identified initial conditions. It is calculated from the maximum daily temperature envelopes and is indicative of the thawing/freezing cycle. Fig. 13B shows the time series of ALT-DOY for the two modelled sites. On average, the propagated uncertainty range in the thaw date is two to four weeks, moving between August and October at the two sites, except for WH’s experiments with dry and average conditions. These configurations showed a large discrepancy at the beginning of simulation up to the beginning of the 2000s and then yielded a similar range of uncertainty.
DZAA (refer to Table 1for definition) is another facet of permafrost that sheds light upon the suitability of the selected soil column depth. It should lie well within the configured depth (e.g . Sapriza-Azuri et al. , 2018), can be used as an auxiliary variable to quantitively assess the simulation (e.g . Elshamy et al. , 2020), or can be used implicitly to identify the presence of permafrost by calculating its corresponding soil temperature (TZAA) (Burke et al. , 2020). In the current study, we assessed the variability in DZAA during the simulation period in response to the initialization. Fig. 13Cpresents DZAA for the whole simulation period (i.e . 1980 - 2016) at the two sites. In general, JMC and WH depict a different response to the spin-up year with average and warm conditions compatible with the simulated ALT and the mean RMSE for the same conditions (seeFig. 11 and Fig. 12 ). The two sites indicate considerable uncertainty to the initial soil storage, with a minor impact of the initial climate condition. In detail, for the JMC site, experiments forced by average and warm climate conditions of the initial year depict an inconsistency of the simulated DZAA, which could reach two-fold (from 6 to 16 m) for tests with intermediate to low initial soil storages (Exps: 11-21). Similar behaviour with higher intensity (four- to-five-fold) is presented at WH for cold and dry initial soil water content cases. The two experiments showed that the selected soil column depth (51.24m) is adequate and suitable to initialize and simulate permafrost for a 38-year period.
Another vital characteristic of permafrost that does not always receive proper attention in the LSMs/ESMs simulations is the permafrost thickness or the depth to the base of permafrost (PB). Modelling PB is only possible in regions with shallow permafrost where it falls within the configured soil depth. PB can be obtained directly from borehole measurements, or from the temperature profile when the thermal probes penetrate all through to the base; derived as the distance between the two (i.e. upper and lower) 0°C isotherms. Alternatively, PB is acquired by extrapolating the temperature envelope using the temperature gradient (i.e.geothermal flux). Besides, PB and MAGTp are among the major factors that describe permafrost’s local and regional conditions and are mainly utilized in engineering design (Wu et al. , 2010), and enhance the simulation/quantification of the global carbon cycle, as being a deep carbon pool (Schuur et al. , 2008). Fig. 13D presents the temporal evolution of the depth to the base of permafrost at the two sites. PB’s general behaviour is comparable to DZAA (Fig. 13C ), with less variability over time, as expected. The simulated permafrost base shows a slight tendency to deepen over time at the two sites under each initializing condition. In other words, the simulations underlined the significant influence of initial soil storage rather than initial climate conditions on the depth to the base of permafrost. Experiments with low initial storage (Exps: 16-21) have the shallowest permafrost base, compared to the intermediate and high values used for initializing soil storage. The simulated PB at the two sites could vary by up to four- to five-fold among all the initialization experiments.
Possible Position of Fig. 13 .