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).