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:
- (subplot A) initial soil
moisture controls spin-up length for temperature regardless of the
spin-up year’s climate. Experiments initialized by a warm climate
showed less dependency on the initial soil moisture.
- (subplot A) the initial year’s climate plays a central role in
determining spin-up length for liquid soil water content. The dry year
requires the most prolonged spin-up to achieve stabilization.
- (subplot A) both the initial year’s climate and initial soil moisture
are comparably important to control the convergence rate of frozen
soil water content.
- (subplot B) the initial year’s climate influences the spin-up needed
for soil temperature. All dry and wet experiments stabilize after less
than 200 cycles, unlike other climate conditions that may need
100-1000 cycles, based on the initial soil moisture.
- (subplot B) the initial climate condition is the primary driver of
liquid water content convergence, similar to the JMC site.
- (subplot B) initial climate and initial soil moisture influence the
required spin-up to attain steady-state frozen soil water content,
with a more substantial influence of the initial climate than the JMC
site.
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:
- Temperature
(subplots A and B): JMC site has a relatively constrained temperature
envelope near the surface (3.5°C), while the biggest variability is
observed at the middle of ODEP (5.5°C) with a diminishing envelope
downward to SDEP (2.5°C). In contrast, the ODEP-SDEP portion of the WH
setup shows the highest variability, varying by 5°C at ODEP and the
range reduces to 1.5°C at SDEP. All JMC experiments were capable of
producing permafrost, which is not the case for the WH setup, as 40%
of its experiments failed to develop permafrost after an exhaustive
spin-up.
- Temperature (subplots A and B):
the impact of initial climate
condition is clear on the simulated temperature at both sites, as
shown by the clustered (well-defined) profiles. Initial soil moisture
partitioning exerts additional influence on simulated temperature as
shown in the scattered profiles of the cold year of JMC site –
moisture partitioning has a limited impact on WH experiments.
- Liquid water content (subplot C): The upper 0.2 m of JMC experiments
(i.e. hemic peat (porosity = 0.93)) has the highest variability
of 0.24 m3 m-3 with a diminishing
envelope downward to ODEP. On the other hand, the whole soil column of
WH site exhibits a wider range of variability for the simulated liquid
soil water profile (0.10-0.15 m3m-3). Still, the surface layer(s) at WH retained
less water, due to the lower soil porosity (0.45). At the bedrock
interface (SDEP), there is a noticeable fluctuation in the simulated
liquid water content for the two setups.
- Liquid water content (subplot C): the initial year’s climate dominates
the amount of stored water in the hemic layer of JMC site, followed by
a negligible impact for the rest of the soil column. Average and warm
years provide the wettest moisture condition for the SDEP portion of
WH. The other climate conditions yield the driest moisture condition
down to SDEP, except the wet year that provides high liquid content
for the upper 2m. The impact of initial soil storage (and the
partitioning) is negligible as experiments with the same climate
condition are distinctively clustered.
- Frozen water content (subplot D): major variability has been recorded
for frozen compared to liquid water content. The two sites show a
similar response across ODEP-SDEP portion (vary by ~
0.40 m3 m-3), accompanied by
significant differences above ODEP. In detail, WH simulations produced
warmer soil temperature (as shown in subplot A and B) down to 1.5-2 m
(ODEP = 0.81 m), while the corresponding depth for JMC case is located
at 0.8 m (within ODEP region of 1.46 m). Such uncertainty in the
frozen water content (entire profile of each site) could propagate
into the simulation and produce divergent results.
- Frozen water content (subplot
D): the two sites did not show any obvious dependency on initial
climate – except warm/average conditions of WH that did not generate
frozen soil, and cold condition in JMC that formed the largest frozen
content at ODEP interface. Experiments with different initial moisture
conditions and the same climate condition are spread within the whole
envelope. This suggests a stronger influence of initial water content
partitioning on the formed ice than the initial climate.Possible Position of Fig. 10 .
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 .