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.