Figure 1. (a) Study area in Eastern Siberia. Black solid line
boxes indicate the imaging areas taken by each satellite. Dashed line
box indicates the area enlarged in (b). Batagay and Verkhoyansk (red
diamonds) are located in the imaging area. (b) Elevation map around
Batagay based on a TanDEM-X DEM (12m mesh). The Batagaika megaslump is
15 km southeast of Batagay. Deformation signals due to the wildfire of
July 2014 were detected in the black dashed area.
Our second objective was to estimate the cumulative spatial distribution
of subsidence, which allows us to estimate the thawed ice volume.
Surface deformation signals over permafrost areas have been interpreted
as being caused by two major processes: (1) irreversible subsidence due
to thawing of ice-rich permafrost or excess ice and (2) seasonally
cyclic subsidence and uplift (Liu et al., 2014, 2015; Molan et al.,
2018). In these previous reports, however, quality interferograms (InSAR
images) were limited in terms of both the temporal coverage and
resolution. This limitation existed because the image acquisition
interval was 46 days at best and the orbit was not well-controlled in
the Japanese Advanced Land Observation Satellite (ALOS) operated from
2006 to 2011 by the Japan Aerospace Exploration Agency (JAXA). For
instance, Liu et al (2015) assumed a simple linear subsidence trend in
their inversion, probably because of the limitation in temporal
coverage. Moreover, the 1.5-year temporal coverage in Molan et al (2018)
would be not long enough to resolve the detailed temporal evolution.
Hence, the total thawed ice volume estimates were uncertain. We also
compared the spatial distribution of subsidence with burn severity and
local landform.
Several studies have reported uplift signals by InSAR over permafrost
areas (Samsonov et al., 2016; Daout et al., 2017; Chen et al., 2018;
Rouyet et al., 2019), but no clear uplift signals have been shown in
previous studies at fire scars as interferometric coherence was lost
during the freezing season in analyzed areas. In contrast, this study
provides the unambiguous detection of upheaval signals in the early
freezing season and confirms the absence of continuing uplift during the
colder season.
Our third objective was, given the clear frost heave signals, to
interpret more physically the observed data. This was because it has
been widely accepted that frost heave is unrelated to volume expansion
of pre-existing pore water into ice, but caused, instead, by ice lens
formation due to the migration of water (Taber, 1929, 1930). However, a
physical understanding of frost heave mechanisms has been established
only during recent decades (e.g., Dash, 1989; Worster and Wettlaufer,
1999; Rempel et al., 2004, Wettlaufer and Worster, 2006; Dash et al.,
2006; Rempel, 2007). Although it appears counter-intuitive, taking a
soil particle inside a unit of ice, there exists an unfrozen (premelted)
water film between the ice and soil even below the bulk-melting
temperature of 0 ℃ (e.g., Dash, 1989; Worster and Wettlaufer, 1999).
Premelted water can be present because of the depression of freezing
temperature by the curved geometry of the soil particle and the
repulsive inter-molecular force between ice and soil particles. Under a
temperature gradient the repulsive thermomolecular pressure on the
colder side is greater than on the warmer side. Hence, the net
thermo-molecular force on the soil particle tends to move it toward the
warmer side, a phenomenon known as thermal regelation (e.g., Worster and
Wettlaufer, 1999; Rempel et al., 2004). Meanwhile, the premelted water
migrates toward lower temperature, where ice lenses will be formed.
These processes are responsible for frost heave and continue as long as
the temperature gradient is maintained, or until significant overburden
pressure is applied (e.g., Dash, 1989; Worster and Wettlaufer, 1999;
Rempel et al., 2004). Although there is still an ongoing debate on the
theory (Peppin and Style, 2013), we applied the simple, physics-based 1D
theory of Rempel et al (2004) to the observed frost heave signal so that
we could physically interpret and explain the observed signals using
reasonable parameters.
2 Study Site
Batagay (67°39′30″ N, 134°38′40″ E) is located on the Yana River, which
is 872 km long and covers a 238,000 km2 basin in a
part of the East Siberian Lowlands in the Sakha Republic (Figure 1). The
elevation ranges between 138 m above sea level at Batagay village and
590 m at Mt. Kirgilyakh on the north-west of Batagaika megaslump (Figure
1b). Our study site was a fire scar located on the western bank of Yana
River, with elevation ~200-400 m (Figure 1b).
The climate is highly continental with a mean annual temperature of
−15.4 °C and mean annual precipitation 170 – 220 mm (Murton et al.,
2017). Meteorological data were sourced from Verkhoyansk, 55 km west of
Batagay. The mean temperature for July and December 2017, respectively,
was 12 °C and −44 °C, while precipitation was 30 mm and 6 mm,
respectively.
We have no in-situ observation data on permafrost conditions and
sedimentology before the fire. However, the burned site is approximately
25 km to the northwest of the Batagaika megaslump (Figure 1); thus, we
refer to the summary provided by Murton et al (2017) as a proxy
for basic information on the burned area and permafrost. The open forest
is dominated by larch with shrubs and lichen moss ground cover. Using
normalized vegetation index by Landsat images we confirmed that the
prefire vegetation at the fire scar was almost the same as that around
the megaslump. Permafrost in the Yana River valley is continuous with
the mean annual ground temperature at the bottom of the active layer,
ranging from −5.5 °C to −8.0 °C, with the active layer thicknesses (ALT)
beneath the forest/moss cover and open sites being 20-40 cm and 40-120
cm, respectively (Murton et al., 2017). In the upslope at Batagaika
megaslump, below the 150 cm thick near-surface sand layer there lies a
20-45 m thick upper ice complex, under which there is a 20-38 m thick
lower sand layer. Below this lies a 3-7 m thick lower ice complex
(Murton et al., 2017). Although the horizontal distribution of this
massive ice complex is yet uncertain, we discuss the possible variations
in the ALT in section 5.2.
The wildfire incident occurred in July 2014 over 36
km2 area, northwest of Batagay (Figure 1). This
wildfire event was evident in the Landsat and MODIS optical images taken
between July 17 and August 2, 2014. While wildfires in northeastern
Siberia are often attributed to human activity (Cherosov et al., 2010),
the cause of the July 2014 wildfire is uncertain. The number of days
with high flammability has noticeably increased over large parts of
Russia, including the Far East (Roshydromet, 2008). For instance, areas
near our study site have experienced even larger wildfires in 2019
(Siberian Times, 2019), as well as a smaller wildfire near the Batagaika
megaslump in 2018.
3 Methods
3.1 InSAR and Data Sets
InSAR has been used as a technique to detect surface displacements (see
Bürgmann et al., 2000; Hanssen, 2001; Simons and Rosen, 2015 for
detailed reviews). InSAR can map surface displacements over the swath
areas with spatial resolution on the order of 10 m or less. InSAR image,
called an interferogram, is derived by taking the differences between
the phase values of SAR images at two acquisition epochs and further
correcting for the known phases contributed from orbital separation
(spatial baseline) and topography. Most SAR satellites have near-polar
orbits, transmit microwave pulses normal to the flight direction and
illuminate the surface of the Earth in ~50-500 km wide
belts depending on satellite type and its observation mode (Figure 1a).
The actual InSAR deformation map indicates the radar line-of-sight (LOS)
changes that are derived by a projection of the 3D surface displacements
onto the LOS direction. Because the incidence angle of the illuminating
microwave is ~30°-40°, LOS changes are most sensitive to
vertical (up-down) displacement followed by east-west displacement and
are least sensitive to north-south displacement because of near-polar
orbit. More specifically, the sensitivity to east-west displacement
changes sign, depending on whether the surface is illuminated from the
east or the west.
Depending on the specific two SAR image pairs and imaged locations, it
is not always possible to quantify surface displacements from
interferograms. As the phase values of an original interferogram are
wrapped into [-π, +π] with 2π ambiguity, they need to be unwrapped
to quantify spatially continuous LOS changes. However, phase unwrapping
becomes impossible when the reflected waves received at the two
acquisitions lack interferometric coherence (i.e., they are uncorrelated
with each other). Lower coherence is caused by long spatial baseline and
temporal changes in the scattering characteristics at the SAR image
resolution cell (temporal decorrelation). For instance, significant
ground cover differences between conditions of deep snow and dry surface
cause temporal decorrelation.
Effects of microwave propagation through non-vacuum medium, ionosphere
and troposphere, on the derived interferometric phase also need to be
considered, as they generate apparent LOS changes that are unrelated to
surface displacements. Moreover, recent studies have also reported the
effect of soil-moisture changes through volume scattering within the
surface soil on the interferometric phase (e.g., De Zan et al., 2014;
Zwieback et al., 2015, 2016).
In this study, we used L-band (23 cm wavelength) HH-polarized SAR images
derived from the PALSAR-2 acquired by the Japanese Advanced Land
Observing Satellite 2 (ALOS2) from 2015 to 2019 together with C-band
(5.6 cm wavelength) VV-polarized SAR images taken during 2017-2019
derived from Sentinel-1 (Figure 2; see also Tables 1 and 2 for details).
The incidence angles at the center of images were 36° and 39° for ALOS2
and Sentinel-1, respectively. In the data sets used, ALOS2 and
Sentinel-1 were illuminating the surface from the west and east,
respectively, and thus the sensitivity to the east-west displacement was
in reverse. To correct for topographic phases, we used TanDEM-X DEM (12m
mesh). Compared to the former ALOS-1/PALSAR-1 InSAR, the ALOS2 orbit is
well controlled, and the spatial baseline is much shorter (Table 1),
which allowed us to ignore DEM errors in the interferograms; the same is
true for Sentinel-1 (Table 2).
Tropospheric delay itself does not depend on the carrier frequency, but
C-band InSAR provides more phase changes because of its shorter
wavelength. In contrast, L-band InSAR phase is more prone to ionospheric
effect, which could be corrected for by range split-spectrum method
(Gomba et al., 2016; Furuya et al., 2017). However, the spatial scale of
ionospheric anomalies was much larger than that of the burned area, and
the ionospheric signals were apparently uncorrelated with the
deformation signal. Thus, we simply took out the long-wavelength phase
trend by fitting a low-order polynomial with clipped InSAR images after
masking out the burned area. We also corrected for topography-correlated
tropospheric errors when they clearly appeared in the InSAR image. These
procedures were somewhat ad-hoc but allowed us to isolate relative
displacements with respect to un-burned areas regarded as reference. It
was also likely, however, that possible long-wavelength permafrost
degradation signals, known as “isotropic thaw subsidence” (Shiklomanov
et al., 2013), were eliminated. Yet, it would be challenging to detect
isotropic thaw subsidence signal only from InSAR data. Hence, we simply
ignored such possible long-wavelength deformation signals.