Possible Position of Fig. 2.
Possible Position of Table 2 .
Possible Position of Fig. 3 .

Model configuration

The JMC and WH models were configured using the approach of Elshamy et al . (2020), where MESH was applied to three permafrost sites along the Mackenzie River valley. Three families of model parameters were identified: soil column configuration, soil texture, and surface canopy. Firstly, for soil layering we employed the scheme of Elshamy et al. (2020) for both sites. This extends to 51.24 m depth and has a fine discretization (9 layers) for the upper 2 meters of the soil, in line with the observed ALT for both sites (Table 3 ). No-heat flux was used as the lower boundary condition of the soil column. The effects of the lower boundary condition were assumed negligible because of its limited impact on the simulated soil temperatures for centennial timescale simulation, as several studies reported (Nicolsky et al. , 2007; Lawrenceet al. , 2008; Hermoso de Mendoza et al. , 2020). Regarding SDEP (depth to the bedrock), we used the gridded bedrock depth dataset by Keshav et al. (2019b) to identify the mean value of SDEP for the two sites: 7.00 m and 10.77 m for JMC and WH, respectively. SDEP is an important parameter that plays a pivotal role in the simulation of permafrost as it alters the thermal regime’s properties (conduction and storage) and the system’s water storage/drainage (see (Elshamy et al. , 2020) for further discussion).
Possible Position of Table 3 .
The second set of parameters are soil texture parameters, which are used in CLASS to parameterize the thermal and hydraulic properties; thus, they strongly influence water and heat storage and movement/conduction. Organic soils are characterized by large heat and moisture capacities that regulate the effects of atmosphere on permafrost around the year (Dobinski, 2011). Therefore, we focused on their representation in our model setups. However, the available soil texture data are insufficient to configure a deep model profile. Soil maps such as Cartographic 1:1000,000 Soil Landscapes of Canada (SLC) v2.2 (Centre for Land and Biological Resources Research, 1996) and its gridded product (Keshavet al., 2019a) only offer data on the spatial characteristics of soil, with no mention of variation with depth. The Global Soil Dataset for Earth system models (GSDE: Wei et al. , 2014) provides gridded texture information for 8 layers but only to a depth of 2.3m. Even though the available borehole logs around (and at) the two sites provide valuable geotechnical data, they lack necessary information on the organic matter content and its thermal and hydraulic properties (the logs provide a qualitative description of the soil components following the USDA classification). Therefore, the most feasible approach was to test different configurations of the soil organic matter, which can be parametrized either as Fully Organic Soil (FOS) or as Mineral Soil with Organic content (MSO). CLASS uses organic matter content within the mineral soil to update only the thermal properties (heat capacity and thermal conductivity), similar to CLM 4.5 (Olesonet al. , 2013). We configured the soil column of JMC using the FOS configuration for the upper 1.46 m (20 cm fibric + 40 cm hemic + 86 cm sapric), and the rest of the soil column as silt loam with high organic content (50% organic content) – full details are given in Elshamyet al. (2020). For the WH site, we configured the upper 0.81 m using the MSO approach (silt loam with 50 - 60% organic content), and the rest of the column as mineral (3.50m sand/silt fine-grained, the rest as silty clay). The MSO configuration was selected for WH site based on the outcomes of benchmarking simulations that yielded better results (i.e. ALT and temperature envelopes) compared to the FOS approach.
The last group of model parameters are those used to parametrize surface canopy conditions. We used the CLASS manual (Verseghy, 2012) to identify the associated parameter values for the two sites, given that JMC is dominated by boreal forest, while evergreen shrubs cover the WH site. Each setup is configured with a single GRU using the single MESH column.

Experimental design

The setups for the two permafrost sites were designed to explore the influence of the uncertain initial model states on the spin-up length (i.e. length required for appropriate model warm-up) and its extended effect on permafrost during a subsequent simulation period. We utilized a single-year spin-up strategy for a maximum of 2000 annual cycles, which is compatible with the available literature on LSM permafrost modelling (Dankers et al. , 2011; Burke et al. , 2013; Elshamy et al. , 2020). As discussed above, the single-year approach is the simplest and most commonly employed method for initializing LSMs (e.g. Rodell et al. , 2005; Nishimuraet al. , 2009; Burke et al. , 2013). However, the resulting state-variables may suffer from an accumulated bias depending on the spin-up year’s climate (Rodell et al. , 2005). Likewise, the other available spin-up techniques (e.g. using a (detrended) sequence of years and long transient simulation) do not resolve this issue (seeSection ‎1 ).
As a single year of forcing is cycled, a complete spin-up would theoretically be achieved if the model states at year m are identical to year m +1. However, achieving a highly precise state equilibrium is not always necessary or feasible, especially for global-scale simulations due to the immense computational cost (Rodell et al. , 2005). We focused on soil temperature profiles and soil moisture (both frozen and liquid) profiles for the stabilization analysis. Previous studies considered either soil temperature (e.g. Burke et al. , 2013; Elshamy et al. , 2020) or total (unpartitioned) soil moisture (e.g. Rodellet al. , 2005; de Goncalves et al. , 2006; Shrestha and Houser, 2010). We considered a tolerance of 0.1°C for temperature states (the same accuracy of permafrost thermal measurements) and 0.01 m3 m-3 for liquid and frozen water states of each soil layer. The stabilization length (i.e. number of spin-up cycles) was defined from when the differences in selected states are less than the identified thresholds, and we used the last timestep of each cycle in the stabilization analysis. There is no consensus among LSMs communities on defining the convergence criteria for adequate model initialization (Yang et al. , 1995). For example, Burke et al. (2013) considered a successful initialization of JULES had been achieved when the variation of soil temperature during spin-up is less than 0.2°C. Employing high thresholds could lead to biased state-variables, while lower thresholds are not feasible for large-scale applications due to their extensive computational cost.
To account for the impact of the initial year’s climate, we selected five climatic conditions based on the total annual precipitation and mean annual air temperature, as per the suggestion by Sapriza-Azuriet al. (2018). We used the WFDEI dataset to identify these, namely wet year (high precipitation), dry year (low precipitation), cold year (low temperature), warm year (high temperature), and an average year (for both precipitation and temperature). Table 4summarizes the climate conditions for the two permafrost sites using a hydrological year (i.e. Oct 1st to Sep 30th).
Similarly, to account for the effects of initial soil moisture content, we considered 21 different uniform cases covering the spectrum of soil water content (water and ice content), as summarized in Fig. 4. These are non-equilibrium states but address the subjectivity of initial soil water content selection/configuration in previous studies. For instance, Rodellet al. (2005) and Shrestha and Houser (2010) defined initial soil moisture as 70% and 10% of saturation for wet and dry conditions, respectively, unlike Cosgrove et al. (2003) and de Goncalveset al. (2006) who quantified these conditions as 100% and 0% saturation, respectively. The two relevant permafrost studies by Sapriza-Azuri et al. (2018) and Elshamy et al. (2020), which utilized the same model as here (MESH/CLASS), did not address this issue. For the dry experiment, with zero saturation for liquid water content, CLASS constrains the residual water content at a value of 4% for mineral soil (MSO configurations). For the FOS configurations, the residual liquid content (or retention capacity) is a function of the organic soil sub-type and varies between 4% and 22%, as mentioned inSection ‎2.1 .
All models were set with the same initial condition for soil temperature, defined as 0°C along the whole profile, except for the bottom temperature, which was extrapolated from the available temperature records following Elshamy et al.(2020), specified as 0.8°C and 0.5°C for the JMC and WH sites, respectively. We assumed uniform initial profiles for soil temperature and moisture contents due to the simplicity of this approach, and to avoid subjectivity in defining a non-uniform profile, especially for temperature. However, the model is shown to rapidly adjust to self-consistent states. Therefore, each site has a total of 105 1-D scenarios {5 climate conditions x 21 soil moisture conditions}, covering distinctive climate and soil moisture conditions. Subsequently, we ran all scenarios for a simulation period of 1979-2016 to assess the impact of uncertainty in initial conditions on various aspects (seeFig. 1 and Table 1 ) of permafrost dynamics. The analysis incorporated quantitative assessment of simulated permafrost in terms of root mean square error (RMSE) for the temperature profiles and ALT error (Bias) over the same period, whenever there were observations.
Possible Position of Table 4.