INTRODUCTION

Permafrost, defined as perennially cryotic ground for at least two consecutive years (Everdingen, 1998), is a critical feature in cold regions that substantially impacts hydrology, energy flux partitioning, plant communities and carbon dynamics. One-half of Canada and one-quarter of the Northern Hemisphere are underlain by permafrost (Zhang et al. , 1999). Permafrost thaw has been observed in recent decades in North America and Eurasia (DeBeer et al. , 2016; Pan et al. , 2016; Meredith et al. , 2020) and is expected to accelerate with global warming (Zhang et al. , 2008a; Lawrence et al. , 2012). CMIP5 results showed a projected areal loss of permafrost in the northern regions of between 30% to 90%, with most of the loss occurring by the end of the 21st century (McGuireet al., 2018). Similarly, Burke et al. (2020) reported projected permafrost degradation of the upper 2m of soil of 10% – 40% per 1.0°C increase in global mean air temperature from the CMIP6 multi-model ensemble. Notably, permafrost contains twice the amount of carbon stored in the atmosphere, the release of which would accelerate the pace of global warming via a positive feedback mechanism (Schuuret al. , 2015; Walvoord and Kurylyk, 2016). Other implications of permafrost thaw include changes to hydrologic connectivity through the formation of vertical and lateral taliks, which in turn affects water fluxes and storage, slope stability and land subsidence, and ecosystem changes including shifts in vegetation structure and streamflow seasonality (Nelson et al. , 2002; Dobinski, 2011; Woo, 2012; Hjort et al. , 2018).
Extensive modelling efforts have focused on permafrost at regional and global scales (Wright et al. , 2003; Lawrence and Slater, 2005; Riseborough et al. , 2008; Zhang et al. , 2013). According to Riseborough et al . (2008), permafrost models can be grouped into three categories. Empirical models are essentially statistical, attribute the occurrence of permafrost to topo-climatic factors (e.g. altitude, air temperature, and slope/aspect) and employ empirically-derived landscape parameters to represent the response of permafrost to climatic and local conditions (Anisimov et al. , 2002; Zhang et al. , 2005). Thus, they are limited to mapping permafrost probabilities, and do not account for thermal inertia (Riseborough et al. , 2008). Equilibrium (or thermal) models use a transfer function between air and ground temperatures to locate the freeze/thaw front. These models,e.g. TTOP (Wright et al. , 2003) and the Geophysical Institute Permafrost Laboratory model (GIPL: Sazonova et al. , 2004), utilize air temperature as their single meteorological input. Such models are applicable only when the system has limited complexity and the transient evolution of permafrost is negligible (Jafarovet al. , 2012; Walvoord and Kurylyk, 2016). Lastly, physically-based models that incorporate heat conduction with phase change are deemed the most accurate method to simulate and project the thermal regime of permafrost over different time periods (Riseboroughet al. , 2008; Jafarov et al. , 2012). The land-surface components of earth system models (ESMs), i.e. LSMs, are one example that couples heat and water transfers across soil-vegetation-atmosphere interfaces.
The current generation of LSMs provides significantly improved simulation of permafrost dynamics. Firstly, the representation of surface insulation by snow cover and soil organic matter has been enhanced, which proved to be a key regulator of atmospheric effects on permafrost thermal and hydraulic regimes (Dobinski, 2011). Improvements include multi-layer snow schemes (e.g. Chadburn et al. , 2015), enhanced frozen soil infiltration algorithms (e.g. Niu and Yang, 2006), explicit organic soil parameterization (e.g. Letts et al., 2000; Lawrence et al ., 2008; Park et al ., 2013), and adding a moss layer as the topsoil layer (e.g. Wu et al. , 2016; Melton et al. , 2019). Secondly, additional permafrost-related physical features/processes have been introduced. Examples include parameterization of hydraulic conductivity to represent the impedance of ice to water movement (e.g. Lawrence et al. , 2011; Wuet al. , 2018), including vegetation dynamics and carbon-pool processes (e.g. Chadburn et al. , 2015; Melton et al. , 2019), and representing lateral taliks (Devoie et al. , 2019) and microtopography (e.g. polygonal regions) (Aas et al. , 2019). Lastly, deeper soils have been used to better capture freeze/thaw cycles (Alexeev et al. , 2007; Nicolsky et al. , 2007). This is significant because the bottom boundary condition strongly dominates soil temperature dynamics at seasonal and longer timescales (Lawrence et al. , 2008) and shallow soil layers are insufficient to resolve the heat storage of the underlying ground. LSM configurations with shallow soil profiles were shown to overestimate the impact of climate change on permafrost (Burn and Nelson, 2006; Burkeet al. , 2020).
Despite these improvements, configuring LSMs to simulate permafrost dynamics is still challenging. The scarcity and uncertainty of permafrost observations limit the representation of spatial heterogeneity over large domains (Chadburn et al. , 2015; Obu et al. , 2019). For instance, the geothermal heat flux is needed to constrain the lower boundary of LSMs; however, no data are available for deep depths and transient conditions (Nishimura et al. , 2009) at large scales. Therefore, most modelling studies either ignore this or assume a constant value spatially and temporally. Moreover, specifying a soil profile of sufficient depth and proper discretization introduces further complications to the modelling. More effort is needed to initialize the deeper profiles due to their long thermal and hydraulic memories (Alexeev et al. , 2007) while the impact of discretization (layering) can be significant and is usually overlooked. In most modelling studies, exponentially increasing soil layer thicknesses with depth are utilized, aiming to balance computational efficiency, numerical stability and model fidelity (Sapriza-Azuri et al. , 2018; Hermoso de Mendoza et al. , 2020). The above factors affect the quality of model initialization and the subsequent simulation. Even in a ‘perfect ’ LSM, improper initialization of state-variables could introduce significant biases in partitioning the surface energy and have a long-lasting effect on model behaviour (Chen and Dudhia, 2001; Rodellet al. , 2005).
Initialization of LSMs can be achieved using observations, despite the scale issues involved, or running model spin-up. In the former, initial conditions are based on realistic field data (e.g. soil moisture, soil temperature). However, gathering extensive datasets over large and remote (e.g. arctic/subarctic) regions and significant depths is challenging (Lamontagne-Hallé et al. , 2020). Alternatively, allowing the LSM to generate its own self-consistent initial conditions through a spin-up technique is commonly used. The spin-up can be performed by repeatedly looping a single year (e.g. CLM: Lawrenceet al. , 2008; NOAH: Shrestha and Houser, 2010; JULES: Dankerset al. , 2011) or a (de-trended) multi-year sequence (e.g.CHANGE: Park et al. , 2013), or through a long transient simulation, of the order of hundreds of years (e.g. CLASS: Sapriza-Azuri et al. , 2018). Looping over a single year is the simplest and most used approach, but may lead to biased estimates of state-variables, depending on forcing (climate) anomalies of the looped year (Rodell et al. , 2005). Similarly, spin-up using a sequence of years is prone to the same problem of bias, in addition to uncertainties associated with de-trending. On the other hand, initial conditions established by running the model for a long transient simulation is no easier. Forcing datasets of sufficiently long periods (100s of years) are rarely available, and can be obtained only through proxy records (e.g. tree rings) or paleoclimatic simulations. Proxy data generally provide mean summer temperature only, are limited spatially, and the reconstructed meteorological variables for LSMs simulation (e.g. precipitation, radiations) suffer from significant uncertainty, as they are typically based on a single variable. Mann et al. (1999) highlighted the non-stationarity in paleo-climatic reconstructions of tree rings during the last millennium, and hence, it is unlikely to find a long enough period of quasi-equilibrium. The simpler approach of single-year spin-up therefore seems to be the most feasible.
In LSM initialization, the relative importance of soil moisture versus temperature memory depends on several factors. Soil moisture memory was shown by Cosgrove et al. (2003) and Rodell et al. (2005) to be larger than thermal memory for shallow soil columns (2m and 3.5m), such that reaching a quasi-equilibrium state for moisture during spin-up requires more time than soil temperature, depending on soil characteristics. However, such conclusions may not be valid for deeper soil columns recommended for LSM simulation of permafrost, due to their larger thermal/hydraulic inertia (Sapriza-Azuri et al. , 2018; Elshamy et al. , 2020; Lamontagne-Hallé et al. , 2020). For example, Elshamy et al. (2020) showed that soil moisture stabilized faster than soil temperature in simulations for the Mackenzie River Basin using a 51.24m soil column. Further, soil moisture data can be helpful. For instance, Walker and Houser (2001) demonstrated added value from assimilating observed surface soil moisture to reduce the spin-up time of a global climate model, noting that different soil properties, e.g. soil texture, hydraulic/thermal conductivity, depth to bedrock, govern the memory of the soil system and hence the required spin-up (Rodell et al. , 2005; Shrestha and Houser, 2010; Elshamy et al. , 2020).
Thus, careful attention to the initialization/specification of soil hydraulic and thermal characteristics is critical for permafrost modelling (Takata, 2002; Lawrence et al. , 2008). Langer et al. (2013) attributed the uncertainties in modelling active layer (i.e. soil depth subjected to seasonal freeze/thaw cycles) dynamics to uncertainties in soil properties and states, especially initial soil water/ice contents. Thermal soil properties depend on moisture content and state (liquid/frozen), especially ice content. Further, the interplay between the climatic conditions of the spin-up year(s) and the initial ice-content of permafrost layers could lead to unrealistic soil thermal and hydraulic properties/states even after reaching an equilibrium state (Rodell et al. , 2005).
Our research addresses two specific objectives: (1) characterizing the impact of initial conditions on the spin-up simulation length required for model warm-up, and (2) quantifying the effect of the uncertainty of model initialization on the simulated permafrost dynamics. Point-scale experiments are configured for two permafrost sites in the Liard River basin. Different combinations of initial soil moisture content (liquid and frozen) and initial year’s climate are used for model initialization. A single-year multi-cycle spin-up strategy is employed for model warm-up, based on the above discussion. The specific contributions of this research are: (1) highlighting the role of initial conditions (climate and soil moisture) for both temperature and moisture and propagating their uncertainty onto a simulation period, and (2) examining various aspects of permafrost dynamics that can be indicative of the quality of the simulation and the parametrization of ground-surface and active-layer. We consider a spectrum of cases for moisture content and its partitioning to liquid and ice, while most previous studies only considered few conditions (Rodell et al. , 2005; Shrestha and Houser, 2010; Burke et al. , 2013). Also, most previous studies selected only a few permafrost characteristics, such as ALT (e.g. Lawrence et al. , 2008; Melton et al. , 2019) or DZAA (e.g. Sapriza-Azuri et al. , 2018; Burke et al. , 2020), to evaluate dynamics (refer toTable 1 for definitions).