Post-wildfire surface deformation at Batagay, Eastern Siberia, detected by L-band and C-band InSAR

Thawing of ice-rich permafrost and subsequent ground subsidence can form characteristic landforms, and the resulting topography they create are collectively called “thermokarst”. The impact of wildfire on thermokarst development remains uncertain. Here we report on the post-wildfire ground deformation associated with the 2014 wildfire near Batagay, Eastern Siberia. We used Interferometric Synthetic Aperture Radar (InSAR) to generate both long-term (1-4 years) and short-term (sub-seasonal to seasonal) deformation maps. Based on two independent satellite-based microwave sensors, we could validate the dominance of vertical displacements and their heterogeneous distributions without relying on in-situ data. The inferred time-series based on L-band ALOS2 InSAR data indicated that the cumulative subsidence at the area of greatest magnitude was greater than 30 cm from October 2015 to June 2019, and that the rate of subsidence slowed in 2018. The burn severity was rather homogeneous, but the cumulative subsidence magnitude was larger on the east-facing slopes where the gullies were also predominantly developed. The correlation suggests that the active layer on the east-facing slopes might have been thinner before the fire. Meanwhile, C-band Sentinel-1 InSAR data with higher temporal resolution showed that the temporal evolution included episodic changes in terms of deformation rate. Moreover, we could unambiguously detect frost heave signals that were enhanced within the burned area during the early freezing season but were absent in the mid-winter. We could reasonably interpret the frost heave signals within a framework of premelting theory instead of assuming a simple freezing and subsequent volume expansion of pre-existing pore water. Post-wildfire surface deformation near Batagay, Eastern Siberia, detected by L-band and Cband InSAR Kazuki Yanagiya and Masato Furuya Department of Natural History Sciences, Graduate School of Science, Hokkaido University. Department of Earth and Planetary Dynamics, Faculty of Science, Hokkaido University. Corresponding author: Kazuki Yanagiya (k.yanagiya@frontier.hokudai.ac.jp) and Masato Furuya (furuya@sci.hokudai.ac.jp)


Introduction
Wildfires in boreal and arctic regions are known to have increased over recent decades in terms of both frequency and areal coverage (e.g., Hu et al., 2010), and have had significant impacts on permafrost degradation (e.g., Zhang et al., 2015;Gibson et al., 2018). Although fires do not directly heat up the subsurface space deeper than 15 cm (Yoshikawa et al., 2003), severe burning decreases surface albedo, and removes vegetation and the surface organic soil layer that previously acted as insulators buffering from changes in air temperature. Subsequent increases in both soil temperature and thickness of the active layer, a near-surface layer that undergoes a seasonal freeze-thaw cycle, have been documented for up to several years after the fire (e.g., Yoshikawa et al., 2003). Meanwhile, in ice-rich permafrost regions, the thawing of permafrost and the melting of massive ice can lead to formation of characteristic landforms such as thaw pits and ponds, and retrogressive thaw slumps. While there are a variety of classifications in terms of morphological and hydrological characteristics , those thaw-related landforms and the topography they create are collectively termed as "thermokarst". However, the role of wildfires in developing thermokarst terrain remains quantitatively uncertain. Moreover, in comparison to the controlled warming experiments in Alaska , wildfires in arctic regions may also be viewed as uncontrolled disturbance experiments that aid in understanding the permafrost degradation processes.
Ice-rich permafrost deposits, known as the yedoma ice complex (yedoma), are widely distributed in the lowland of Alaska and Eastern Siberia . The greatest subsidence within the 2007 Anaktuvuk River tundra fire scar was identified in the yedoma upland by LiDAR . Yedoma is a unique permafrost deposit in terms of its extraordinarily high volume of ice (50-90 %) and organic-rich sediments. While the organic carbon trapped in permafrost regions is estimated to be twice that in the current atmosphere, permafrost thawing and related thermokarst processes may release the carbon as greenhouse gasses (CO 2 and CH 4 ) via microbial breakdown, which may further promote global warming . Thus, in order to estimate the volume of greenhouse gasses released, it is important to evaluate the volume of thawed ice associated with thermokarst processes in yedoma-rich areas.
Near the village of Batagay, Sakha Republic, Eastern Siberia (Figure 1), there exists the Batagaika megaslump, known as the world's largest retrogressive thaw slump, exposing roughly 50-90 m thick yedoma deposits on the north-east facing slope (e.g., . Thaw slumps are characterized by a steep headwall surrounding a slump floor and develop as a result of rapid permafrost thawing. The Batagaika megaslump was initiated at the end of the 1970s by deforestation but still appears to be growing . Considering this feature, it is worth considering whether new disturbances in the proximity will result in the formation of similar landforms. A wildfire incident occurred in July 2014 near Batagay, which, like deforestation, will change the ground thermal regime. Therefore, it is important to examine whether future catastrophic thermokarst development could be similarly initiated at the fire scar, whose area is much larger than the Batagaika megaslump (Fig 1b).
The first objective of this study was to assess the effectiveness of satellite Interferometric Synthetic Aperture Radar (InSAR) in detecting surface deformation signals due to wildfire-induced thermokarst over different temporal scales. InSAR has been used to detect long-term and seasonal displacements over several thawrelated landforms in permafrost areas (e.g., Iwahana et al., 2016;Antonova et al., 2018;. Although subsidence signals as a result of thermokarst associated with Alaskan wildfires have been detected using InSAR Iwahana et al., , 2016b, no such studies have been conducted on Siberian fires, to our knowledge. Also, all previous InSAR-based post-wildfire deformation mapping has been performed over relatively flat terrains, but no reports over hillslopes have been shown. Moreover, in contrast to previous studies, we employed two independent SAR imageries with distinct carrier frequencies and polarizations, L-band (1.2 GHz) HH-and C-band (5.4 GHz) VV-polarized microwave. Because the imaging geometries were different and had different sensitivities to the 3D displacement vector, we could not only take advantage of the performance of each sensor in mapping deformation signals but could also cross-validate the measurements by two InSAR data sets. 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 . 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,  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  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 , 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 . However, a physical understanding of frost heave mechanisms has been established only during recent decades (e.g., . 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., . Premelted water can be present because of the depression of freezing temperature by the curved geometry of the soil particle and the repulsive intermolecular 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., . 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., . Although there is still an ongoing debate on the theory , we applied the simple, physics-based 1D theory of  to the observed frost heave signal so that we could physically interpret and explain the observed signals using reasonable parameters.

Study Site
Batagay (67°39'30" N, 134deg38'40" E) is located on the Yana River, which is 872 km long and covers a 238,000 km 2 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 degC and mean annual precipita-4 tion 170 -220 mm . Meteorological data were sourced from Verkhoyansk, 55 km west of Batagay. The mean temperature for July and December 2017, respectively, was 12 degC and -44 degC, 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  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 degC to -8.0 degC, with the active layer thicknesses (ALT) beneath the forest/moss cover and open sites being 20-40 cm and 40-120 cm, respectively . 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 . 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 km 2 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 , 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 . 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.

InSAR and Data Sets
InSAR has been used as a technique to detect surface displacements (see Burgmann 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˜30deg-40deg, 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 [-π, +π] ωιτη 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 5 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., 2015Zwieback et al., , 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°f or 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 . 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" , 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. While ALOS2 has only imaged the area since 2014 its data acquisition interval is much longer than that of Sentinel-1 (Figure 2). Previous studies demonstrate that it is not possible to infer the total subsidence using pre-and post-wildfire SAR images, as the drastic changes in land cover cause low interferometric coherence . Additionally, JAXA changed the carrier frequency of PALSAR-2 in June 2015 ( Figure 2). Hence, monitoring long-term deformation using ALOS2 InSAR is possible only since October 2015. Conversely, frequent data acquisition in Sentinel-1 started only in 2017. Thus, we first performed an inter-comparison between ALOS2 and Sentinel-1 InSAR, focusing on the seasonal changes in 2017. We stacked Sentinel-1 interferograms to set the temporal coverages to nearly identical with those of ALOS2 ( Figure 2). Stacking was necessary because we failed to derive long-term Sentinel-1 interferograms, as the burned areas quickly lost coherence. During the temporal interval of ALOS2 images Sentinel-1 had more cycles. Therefore, the number of Sentinel-1 stacks varied from three to eight (Figure 2).
Although L-band SAR is known to have better interferometric coherence than C-band SAR (e.g.,  our study indicated that Sentinel-1 could maintain a comparable interferometric coherence with L-band ALOS2 even during winter season. This is likely due to the short acquisition period of 12 days as well as the somewhat drier snow in the area in winter that allows microwave to reach the ground. A dry snow cover of depth less than 1 m is undetectable to microwave radiation, whereas over wet snow surface scattering dominates . The frequent data acquisition of Sentinel-1 since 2017 allowed us to examine detailed seasonal changes in surface deformation ( Figure 2). Some Sentinel-1 InSAR pairs in early summer, however, did not show good coherence, possibly due to snow wetness.
In order to infer long-term temporal changes and cumulative displacements, we performed SBAS (Small Baseline Subset)-type time-series analysis Schmidt and Bürgmann, 2003), using 50 high-quality ALOS2 interferograms that included one-year-as well as short-term interferograms ( Figure  3). We could estimate the average LOS-change rates between each acquisition epoch without assuming any temporal change models. In contrast to the original SBAS approach, we did not estimate DEM errors because the well-controlled orbit, as well as the precise TanDEM-X DEM, have no sensitivities to those errors.
In order to estimate the errors of the derived time series, we assumed each original SAR scene contained 0.2 cm errors, and made InSAR data covariance matrix, following the method of . The errors are relatively smaller than those in previous studies of SBAS analysis (e.g. 0.4 cm in 0.75 cm in Biggs et al., 2007), because, as noted earlier, we took out the long-wavelength phase trend form each InSAR image, and the analysis area is smaller (12× 12 km) than previous studies.

Multispectral remote sensing of burn severity
Normalized burn ratio (NBR) is a useful multispectral remote sensing index to assess the impact of wildfire on vegetation. Vegetation reflects more strongly in the near-infrared (NIR) than in the shortwave infrared (SWIR) region, while a fire scar reflects more strongly in the SWIR. Utilizing this property, NBR is defined as NBR=(NIR-SWIR)/(NIR+SWIR). The difference NBR (dNBR), i.e., the difference between prefire NBR and postfire NBR, indicates burn severity (Key and Benson, 2006;Miller and Thode, 2007). Generally, when dNBR is greater than 0.66 the fire is regarded as "highly severe". We computed dNBR for the 2014 fire using Landsat 8, Band 5 (850-880 nm) and Band 7 (2110-2290 nm) images for near-infrared and shortwave-infrared, respectively, to associate the inferred subsidence distribution with burn severity.

One dimensional frost-heave theory based on premelting dynamics
We used the one-dimensional frost-heave theory as a tool to interpret the observed uplift signals. Inspired 8 by one-way frost heave experiments ,  and  derived a steady-state heave rate V l of an ice lens, considering the force balance among thermo-molecular force F T , hydrodynamic force F µ , and overburden force F O (pressure P 0 ). Here, we assumed a constant heave rate V l , which may not necessarily reflect the actual observations shown below as well as in . However, this assumption simplified the theory, and we assumed that the observed heave rate did not change drastically over time.  proposed a non-dimensional heave rate v l of an ice lens as a function of its boundary position ξ l given: where µ, k 0 , and ρ are the viscosity of water, the permeability of ice-free soil, density of water, respectively. The quantityG ≡ L Tm ∇T has the same dimension as gravity and indicates thermo-molecular force when multiplied by the mass of displaced ice; L is the latent heat of fusion and T m is the bulk melting temperature. The first and second term in the bracketed numerator are proportional to F T andF O , respectively, while the bracketed denominator is proportional to F µ . The integral is performed alongξ ≡ z z f , where z f is the position above (below) where ice saturation S s becomes non-zero (zero);z h indicates the position where hydrostatic pressure is achieved, and φ is the porosity of soil. The normalized overburden pressure and permeability are defined asp 0 ≡ P0 ρΓz f andk ≡ k k0 ≥ 1, respectively.

Results
We performed an inter-comparison of ALOS2/Sentinel-1 interferograms, focusing on the seasonal changes in surface deformation. We then showed short-term deformation derived by Sentinel-1 and long-term deformation derived by time-series analysis of ALOS-2. Subsequently, we estimated the total volume of thawed excess ice. Although both satellite images covered the Batagaika megaslump we did not observe clear LOS changes as detected at the fire scar, which could be due to the lack of spatial resolution of the InSAR images.  (Bottom) Sentinel-1 stacked interferograms during the five periods, (f)-(j), derived so that the temporal coverage could nearly match those from (a) to (e); all the interferograms are overlaid on shaded relief maps. Warm and cold colors indicate LOS changes away from and toward the satellite, respectively.

9
We compared the ALOS2 and stacked Sentinel-1 interferograms for five periods ( Figure 4) and assessed their differences ( Figure 5). Despite differences in look directions both ALOS2 and Sentinel-1 indicated extensions in the LOS during periods (a, f) from the middle of June to the end of July and (b, g) from the end of July to the early October. Also, their deformation areas and amplitude were mostly consistent, suggesting that LOS changes were largely due to summer subsidence (see section 5.1 below for details). In terms of the spatial distribution of deformation signals, we noticed that the LOS changes over higher-elevation areas such as ridge and peak were insignificant, whereas the boundaries between the burned and unburned areas were clear. The north-western area, however, showed few LOS changes (see section 5.2 for the relationship between LOS changes and burn severity). During the period (c, h) from early October to early December both ALOS2 and Sentinel-1 indicated shortening in the LOS by an approximate 5 cm maximum, and the deformation areas and amplitude were quite similar. This observation presumably indicated frost heave in the early freezing period. In view of the previous three periods, both subsiding and uplifting areas were nearly the same. The following period (d, i) from early December to the middle of March also included the winter season with much colder air temperatures, but we did not observe any significant deformation signals, indicating that frost-heave virtually stopped in early December. While the good interferometric coherence during mid-winter was an unexpected result, we speculate that it could have been due to drier, lower amounts of snowfall.
In the periods (e, j) from the middle of March to early June, both ALOS2 and Sentinel-1 suffered from decorrelation, and we could not identify clear deformation signals. However, in light of Figure 6 below, each of the Sentinel-1 interferograms had overall good coherence with the exception of the data acquired in the middle of May. These observations suggested that the decorrelation may be attributable to the rapid changes on the ground surface during the initiation of the thawing season when the air temperature rises above the freezing point and the active layer begins to thaw.  (Figure 4e and 4j), we could not estimate differences due to coherence loss. Figure 5 shows the differences between ALOS2 and Sentinel-1 InSAR data with nearly identical periods, which may help in cross-validating the measurements and understanding the actual deformation processes. The estimated differences and their 2σ scatter were 0.5±1.2cm (Fig 5a),0.7±2.3cm (Fig 5b),0.3±1.3cm (Fig  5c), and 0.6±0.3cm (Fig 5d), with mean of 0.5±1.5 cm. The differences and their variances were variable over time but apparently indicated some systematic trends. For instance, over the east-facing slopes, the differences were almost always positive (This is discussed more comprehensively in section 5.1).
The Sentinel-1 interferograms for 2017 demonstrate that the progress of deformation was not at a constant rate ( Figure 6). The most rapid deformation took place in June (periods 1 and 2) with no substantial deformation in July (period 3) and started to subside again in August (periods 4-6). We found that the subsidence occurred sporadically over time and space and that the burned area did not uniformly subside. For periods 4, 5 and 9 we were unable to perform phase unwrapping at specific locations near the ridge and the boundaries between the burned and unburned areas. These unwrapping errors were responsible for the localized, large differences observed in Figure 5b. We confirmed the presence of low coherence bands along the unwrapping errors, which may suggest large phase jumps due to large displacements during the 12 days; enigmatically, no such line-shaped low coherence was detected in the long-term ALOS2 interferograms. Moreover, Figure 6 demonstrates that the frost heave started in late September, which was missed in the periods (b) and (g) of Figure 4, and that the absence of any deformation signals lasted from early December to May of the following year. We will physically interpret the absence of deformation signals during the coldest season in section 5.3.  Table 2.  Table 1; imaging was performed by ascending, right-looking orbit. Warm and cold colors indicate LOS changes away from and toward the satellite, respectively. Black dashed line indicates the boundary between the burned and unburned area confirmed with Landsat optical images.

Long-term deformation inferred from time-series analysis of ALOS2 interferograms
Figures 7a-7h show ALOS2 interferograms, each of which covers nearly one-year after October 2015 with some overlaps in its temporal coverages. Figure 7a, derived at the earliest period after the fire, indicates the maximum one-year subsidence to be as much as 10 cm or more.
If the amplitude and timing of seasonal subsidence/uplift cycle are invariable over time, a one-year interferogram will tell us only the irreversible displacements regardless of the acquisition times of master/slave images, which corresponds to the "pure ice" model in . Figure 7 sequentially shows the periods from October 2015 to June 2019 and indicates that the yearly subsidence rate slowed down. However, the variations of the one-year LOS changes in Figures 7 suggest that the actual deformation processes were more complex. Figure 8a shows the cumulative LOS changes from October 2015 to June 2019 derived from SBAS-type timeseries analysis, and that the maximum LOS extension reached as much as 25 cm; the 2σ errors for Figure 8a were ±1.5cm. Considering that the LOS changes during the first year after the 2014 fire were not included, the total LOS changes were presumably much greater than 25 cm, which meant that the subsidence was greater than 30 cm on account of the 36°incidence angle. As mentioned earlier, however, the higher-elevation areas such as the ridge did not undergo significant deformation, which probably would have been the case even during the first year after the fire. In addition to the high elevation areas, we realized clear contrasts in the LOS changes between the east-and the west-facing slopes near the northwestern area and the central north-south trending ridge; this spatial heterogeneity could also be recognized in Sentinel-1 ( Figure 4). Their possible mechanisms comparing the burn severity ( Figure 8b) and local landform ( Figure 8c) are discussed in section 5.2.
We show the estimated time-series data at four representative sites (Figures 9a-9d), whose locations are indicated in Figure 8a. The sites (a) and (b) underwent nearly the same cumulative LOS changes by roughly 20 cm but were located at different slopes that are 4.3 km apart. On the other hand, the cumulative LOS changes at the site (d) were relatively small (approximately 10 cm). The site (c) located in the ridge did not show either significant seasonal or long-term deformation.
Time series data in Figures 9a and 9b clearly indicate that the largest subsidence took place from 2015 and 2016. We believe, however, that the most significant subsidence probably occurred only during the thaw season in 2016, as we have observed earlier, that no deformation occurred from December to March. Thus, the actual subsidence rate from October 2015 to July 2016 should have been more complicated than that expected from the linear trend in Figures 9a and 9b. The error bars in Figures 9a-9d indicated an estimated standard deviation with 2σ and attained ±1.5cm in the last epoch.

Estimating the total volume of thawed excess ice
Post-wildfire deformation over a permafrost area presumably consists of two contributions: (1) irreversible subsidence due to melting of ice-rich permafrost below the active layer, and (2) seasonally cyclic subsidence and uplift due to freeze-thaw of the active layer . In order to separate the two processes from the observed deformation data,  used independent groundmeasured ALT data to predict the ALT contribution to total subsidence. Ground-measured pre-fire ALT data were not available at this study site. Given the temporal evolution of post-wildfire deformation data (Figures 9a-9d), however, we may regard the cumulative deformation in Figure 8a as being due to irreversible subsidence during the period between October 2015 and June 2019, and estimate the total thawed volume as 3.56 ± 2.24 × 10 6 m 3 ; the error bar is based on the root mean square of the no-deformation signals outside the burned area, which is multiplied by the burned area. However, in view of the temporal evolution in Figure 9, we could speculate that a much larger deformation was also taking place immediately after the 2014 fire until October 2015, during which, unfortunately, no deformation data are available. Thus, this estimate should be viewed as a lower estimate, with the actual volume of thawed permafrost possibly being much greater.
Nevertheless, despite its much smaller area size (Figure 1b), the thawed volume at the Batagaika megaslump is greater than 2.5× 10 7 m 3 , an order-of-magnitude larger than our estimate above. Moreover, the thaw-subsidence rate at the fire scar is slowing down (Figure 9). We discuss the possibility of the another megaslump emergence at the fire scar in section 5.2. 5 Discussion 5.1 Similarities and differences in the ALOS-2 and Sentinel-1 interferograms: implication for insignificant slope-parallel sliding Taking into account the imaging geometries of the ALOS2 and Sentinel-1, we could comprehensively interpret the differences in Figure 5 and also infer the actual deformation processes. The weights multiplied to the 3D-displacements, (U ew , U ns , U ud ), to compute LOS changes were +0.573, +0.132, and -0.809 for ALOS2 and -0.583, +0.236, and -0.777 for Sentinel-1, respectively; eastward, northward and upward displacements were taken to be positive. Assuming the LOS changes of the two sensors are identical (which is roughly the case in Figure 4), and no north-south displacement U ns , the constraint on the east-west and up-down displacements can be derived as U ew :U ud = 0.032:1.156. The assumption of zero U ns might appear unrealistic but can be reasonable over the east-and west-facing slopes, which incidentally cover a broad area of the fire scar. As this constraint indicates the dominance of vertical displacement, we can infer that slope-parallel sliding did not take place over the east-and west-facing slopes.
In the thawing season when the vertical displacement is downward (negative), the previous constraint on the two displacements also indicates that the east-west displacement should always be westward (negative), regardless of the slope. As this is physically implausible, we may assume that east-west displacements were virtually zero over both the east-and west-facing slopes. We can thus infer a pure vertical subsidence without any east-west displacements during the thaw season. Hence, the differences between ALOS2 and Sentinel-1 in the thawing seasons (Figures 5a and 5b) will be simply equal to -0.032U ud . Therefore, we can expect systematically positive differences in the thawing season, regardless of the east-and west-facing slopes, which appear consistent with observations (Figures 5a and 5b). Quantitatively, however, the mean differences of 0.5-0.7 cm are too large to be attributable to the geometric difference alone, on account of the subsidence by as much as 5 cm or more. Here, we hypothesize the possible impact of soil-moisture changes, which can reach˜10 % of the carrier wavelength (Zwieback et al., 2015(Zwieback et al., , 2016. As changes in soil moisture generate larger apparent LOS changes in L-band than in C-band InSAR, the observed differences can be likely. In contrast to thaw subsidence frost-heave is more likely to occur towards the slope normal direction. Assuming the magnitude of slope-normal uplift, U f , over a slope with gradientθ , the differences between ALOS2 and Sentinel-1 would be U f (1.156 sinθ -0.032cosθ) assuming zero U ns . We estimated |θ |=1.58°, which corresponds to 55 m height difference over 2 km horizontal distance and was fairly consistent with the slope of the studied area. Meanwhile, the differences can also be considered 1.156U ew -0.032U ud , which indicates additional positive and negative effects on the east-and west-facing slope, respectively. Indeed, Figure 5c appears to depict clearer contrasts in sign on the east-and west-facing slopes. Moreover, the impact of changes in soil moisture are likely much smaller in the colder season than in the thaw season, which may explain the smaller differences in freezing seasons.
As we derived the differences over the no-deformation season (Figure 5d), we can attribute them to the at-mospheric effect on ALOS2 and Sentinel-1 interferogram (Figures 4d and 4i). The overall positive differences are likely because the spatial scale of atmospheric delay was greater than the fire scar area.
Previous reports of thermokarst subsidence after fire have focused on relatively flat areas as those at the 2002 tundra fire in the central Seward Peninsula, Alaska (Iwahana et al., 2016b), the 2007 Anaktuvuk River tundra fire , and the 2009 Big Creek Fire in the Alaskan Yukon River basin . As such, in addition to the broad subsidence detected by InSAR, polygonal patterns associated with ice wedge degradation became clearly visible 4-7 years after the fire by high-resolution optical and LiDAR remote sensing Iwahana et al., 2016b). At the studied hillslopes, in contrast, no such polygonal patterns are likely to be detected. Nonetheless, the dominance of vertical displacements with little slope-parallel sliding indicate that rapid active-layer detachment sliding (ALDS) events were insignificant. In contrast, many ALDS events triggered by fire have been mapped at Mackenzie Valley, Canada, whose length could sometimes reach hundreds of meters . If ALDS event with such length occurred, we could have observed significant loss of interferometric coherence. It is possible, however, that local ALDS events occurred but were undetected because of the coarse resolution (˜10 m) of InSAR images. Because the subsidence was caused by thawing of ice-rich permafrost, meltwater should have been supplied at the base of active layer. Considering the mechanisms of ALDS , porewater pressure increase might have been not enough to reduce the effective overburden stress and to initiate significant slope-parallel sliding. This is possibly because the meltwater could have drained through the gullies. However, in view of the significant uplift signals over the burned area even years after the fire, the meltwater is still likely to be undrained on the slope. If there were further enough water input by, for instance, warmer days and/or heavy rain, significant ALDS events may take place in the future.

What controls the heterogeneous distribution of subsidence magnitude? Possible emergence of another megaslump
The cumulative subsidence magnitude was spatially variable but showed some systematic changes. In addition to the ridges and peaks the west-facing slopes showed significantly smaller subsidence than the east-facing slopes (Figures 4 and 8a). To interpret the spatially heterogeneous subsidence, we associated burn severity and local landform with the cumulative subsidence ( Figure 8). In light of the inferred dNBR (Figure 8b), which ranged from 0.2 to 0.4, the burn severity was moderate rather than high. Also, the burn severities were spatially less heterogeneous than those of cumulative subsidence and local landform. In fact, we could even identify deformation-free areas having even higher burn severity. Thus, although the fire undoubtedly initiated the subsidence, the burn severity did not control the subsequent cumulative magnitude.
Notably, however, gullies were clearly more developed on the east-facing slopes than on the west-facing slopes (Figure 8c), which were confirmed to be present at least back in 1991 by Landsat image. Considering the striking correlation between the development of gullies and the larger subsidence, there is high likelihood of a causal relationship between them. Similar dependence on the slope aspect was reported by Lacelle et al (2010,2015), who found that hillslope thaw slumps in the Richardson Mountains-Peel Plateau, northwest Canada, predominantly developed on the east-facing slope.  interpreted that the active layer on the east-facing slope might be thinner because of lower amount of insolation than on the south-and west-facing slopes, which promoted a triggering mechanisms of thaw slumps because the ice-rich permafrost was closer to the surface. Although the broadly subsiding areas are not so-called thaw slumps, thinner active layers on the east-facing slope are likely and can consistently explain both the larger subsidence and the rich development of gullies. This hypothesis can be tested either by examining the surface deformations at the 2018 and 2019 fire scars and other fire scars across Siberia and other boreal regions or by performing field-based thaw-depths measurement.
The recent slowdown of the subsidence rate ( Figure 9) may suggest that the 2014 fire scar could stabilize in the near future. However, although it depends on how quickly the vegetation is recovered, we do not preclude the possible emergence of another megaslump particularly on the east-facing slopes. In order to initiate thaw slumps, ice-rich permafrost needs to be exposed at the surface , at which the initial headwall and slump floor are formed. In contrast to the thaw slumps near shorelines, coastlines and riverbanks (e.g., Kokelj et al., 2009), no mechanical erosions by waves and currents are effective on hillslopes like the studied area. For the development of retrogressive thaw slumps (RTS) on hillslopes,  suggested that ALDS triggered by meteorological events could remove the overlying active layer and expose the ice-rich permafrost. Although no large-scale ALDS events were detected during the studied period, they might take place as discussed in the previous section. Moreover, Figure 8a indicates that the subsidence magnitude becomes larger toward upslope, and there are clear boundaries between the subsiding and non-subsiding portions, where an initial headwall for RTS could be exposed. Once an initial headwall has formed, subsequent retreat rate is rapid on the order of several meters per year . Thus, in order to monitor the early formation process of RTS in more detail, long-term radar remote sensing with higher spatial and temporal resolution would be necessary and promising. Figure 10. Non-dimensional heave rate profiles of an ice lens as a function of its boundary position, based on the analytical model by . Five cases of non-dimensional overburden pressure p 0 and porosity φ are shown.

Interpretation of frost-heave signals
In order to interpret the detected frost heave signal on the basis of the theory by , we first examine the sensitivity of the heave rate on the normalized overburden pressure p 0 and porosity φ. Figure  10 shows five cases of non-dimensional heave rate profiles as a function of the ice lens boundary positionξ l , indicating that the maximum heave rate is mainly controlled by the normalized overburden pressure p 0 and is somewhat insensitive to the porosity φ. Details of the heave rate profiles will depend on the assumed models of permeability and ice saturation, but the qualitative characteristics are not altered . There exist two positions that give the same heave rate, but only the branch with smaller ξ l is stable .
We can attribute the clear contrast in the frost heave signals inside and outside the burned area to the differences in the normalized overburden pressure p 0 . Because the mechanical overburden pressure P 0 will not significantly differ from the inside to the outside of the burned area, the larger frost heave rate in the burned area would be caused by larger temperature gradient G and/or deeper frozen depth z f . Owing to the removal of vegetations and surface organic layers over the burned area, the larger temperature gradientG than that of the unburned area is likely more marked in the early freezing season and may generate a greater thermomolecular force that will effectively reduce the normalized overburden pressure. We may also interpret the absence of frost heave signals in mid-winter as due, probably, to the smaller temperature gradient G than that in late fall/early winter; if frost heave were controlled by temperature instead of temperature gradient, we would expect even more significant signals during the much colder part of the season. The deeper frozen depthz f is also likely due to the loss of surface vegetation and should supply more water for frost heave.
From the end of September to the middle of November 2017, Figure 6 shows LOS changes by approximately 1.5 cm over 12 days toward the satellite that corresponds to an approximate 1.9 cm uplift. Assuming a constant-rate frost heave, this corresponds to a heave rate of 1.8× 10 -8 (m/s). The most critical parameter controlling heave rate is the permeability for ice-free soil k 0 , which can vary by orders-of-magnitude, while other parameters are well-constrained. We may fit our observed heave rate with the ice-free permeability, k 0˜1 0 -17 (m 2 ), which is a likely value in view of the three cases in .
Here we comment on the modeling of uplift signals as caused by in situ freezing of pore water into ice . The in situ freezing model is simple, and can explain the timing, duration, and magnitude of uplift signals, if one assumes such pore water in the active layer. However, because the Stefan function approach in Hu et al (2018) is essentially controlled by atmospheric (or ground) temperature that is rather homogeneous over this spatial scale, it is difficult to account for the observed heterogeneous distribution of uplift signals. The distribution of uplift signals was closely correlated with that of subsiding signals, which led us to interpret that the permafrost thaw and its incomplete drainage could become a water reservoir for ice lens formation and frost-heave. The frozen pore ice within the soil and the ice lens formed by water migration are totally different in terms of their formation mechanisms and subsequent forms of ice. From a geomorphological perspective, the presence of ice lenses will play a role in reducing the strength of soil and potentially initiating ALDS, because porewater pressure will increase at the front of thawing, whereas pore ice within the soil would simply stay as pore water with little impact on the landform.
We also recognize, however, that the microphysics-based theory adopted in this study is developed in 1-D geometry and is based on the assumption of "frozen fringe", a region where liquid freezes into ice through the pores of soil. Some laboratory experiments did not support the presence of frozen fringe (e.g., , and the "fringe free" frost heave theory has also been proposed; see  for review. In addition to the controlled lab experiments and theoretical developments, more detailed observations of natural frost heave signals are becoming possible and might help better understand the physics of frost heave and its geomorphological consequences.

Conclusions
We used L-band and C-band InSAR to detect post-wildfire ground deformation at Batagay in Sakha Republic, showing not only subsidence signal during the thawing season, but also uplift during the early freezing season and virtually no deformation in midwinter without loss of coherence. Time series analysis allowed us to estimate cumulative displacements and their temporal evolution, as quality interferograms could be obtained even in the winter season. We found that the thawing of permafrost in the burned area lasted three years after the fire, but apparently slowed down after five years. During the studied period, no significant slope-parallel sliding was detected, and the post-wildfire deformation was mostly subsidence. Despite the rather homogeneous burn severity, the cumulative subsidence magnitude was larger on the east-facing slopes and showed a clear correlation with the development of gullies, suggesting that the east-facing active layers might have been originally thinner. Short-term interferograms (2017Short-term interferograms ( -2018 indicated that the subsidence and uplift was clearly enhanced compared with the unburned site. We have thus interpreted the frost heave signals within a framework of premelting dynamics. Post-wildfire areas are a focus of permafrost degradation in the Arctic region. This study is supported by Researcher's Community Support Projects of Japan Arctic Research Network Center in 2016-2019, andby KAKENHI (19K03982). PALSAR2 level 1.1 data are provided by the PALSAR Interferometry Consortium to Study our Evolving Land Surface (PIXEL) and the ALOS2 RA6 project (3021) under cooperative research contracts with the JAXA. Sentinel-1 SLC data are freely available.
TanDEM-X DEM copyrighted by DLR and were provided under TSX proposal DEM GLAC1864. Climate data at Verkhoyansk, Russia, are available from ClimatView site; http://ds.data.jma.go.jp/tcc/tcc/products/climate/climatview/outline.html. We thank Go Iwahana for discussing our preliminary results. We also acknowledge Lin Liu, two anonymous reviewers and the editors, Joel B. Sankey and Amy East, for their extensive and constructive comments, which were helpful in improving the original manuscript.  Spatial heterogeneities of the subsidence magnitude were clearly correlated to the gully development, whereas the burn severity was rather homogeneous.
 Detection of enhanced uplift signals at the fire scar and its interpretation based on a physics-based frost heave theory.

Abstract
Thawing of ice-rich permafrost and subsequent ground subsidence can form characteristic landforms, and the resulting topography they create are collectively called "thermokarst". The impact of wildfire on thermokarst development remains uncertain. Here we report on the post-wildfire ground deformation associated with the 2014 wildfire near Batagay, Eastern Siberia. We used Interferometric Synthetic Aperture Radar (InSAR) to generate both longterm (1-4 years) and short-term (sub-seasonal to seasonal) deformation maps. Based on two independent satellite-based microwave sensors, we could validate the dominance of vertical displacements and their heterogeneous distributions without relying on in-situ data. The inferred time-series based on L-band ALOS2 InSAR data indicated that the cumulative subsidence at the area of greatest magnitude was greater than 30 cm from October 2015 to June 2019, and that the rate of subsidence slowed in 2018. The burn severity was rather homogeneous, but the cumulative subsidence magnitude was larger on the east-facing slopes where the gullies were also predominantly developed. The correlation suggests that the active layer on the east-facing slopes might have been thinner before the fire. Meanwhile, C-band Sentinel-1 InSAR data with higher temporal resolution showed that the temporal evolution included episodic changes in terms of deformation rate. Moreover, we could unambiguously detect frost heave signals that were enhanced within the burned area during the early freezing season but were absent in the mid-winter. We could reasonably interpret the frost heave signals within a framework of premelting theory instead of assuming a simple freezing and subsequent volume expansion of pre-existing pore water.

Plain Language Summary
Wildfires in arctic regions not only show an immediate impact on nearby residents but also long-lasting effects on both regional ecosystems and landforms of the burned area via permafrost degradation and subsequent surface deformation. However, the observations of post-wildfire ground deformations have been limited. Using satellite-based imaging technique called Interferometric Synthetic Aperture Radar (InSAR), we detected the detailed spatial-temporal evolution of post-wildfire surface deformation in Eastern Siberia, which helps in understanding permafrost degradation processes over remote areas. Post-wildfire areas are likely to be focal points of permafrost degradation in the Arctic that can last many years.

Introduction
Wildfires in boreal and arctic regions are known to have increased over recent decades in terms of both frequency and areal coverage (e.g., Hu et al., 2010), and have had significant impacts on permafrost degradation (e.g., Zhang et al., 2015;Gibson et al., 2018). Although fires do not directly heat up the subsurface space deeper than 15 cm (Yoshikawa et al., 2003), severe burning decreases surface albedo, and removes vegetation and the surface organic soil layer that previously acted as insulators buffering from changes in air temperature. Subsequent increases in both soil temperature and thickness of the active layer, a near-surface layer that undergoes a seasonal freeze-thaw cycle, have been documented for up to several years after the fire (e.g., Yoshikawa et al., 2003). Meanwhile, in ice-rich permafrost regions, the thawing of permafrost and the melting of massive ice can lead to formation of characteristic landforms such as thaw pits and ponds, and retrogressive thaw slumps. While there are a variety of classifications in terms of morphological and hydrological characteristics , those thaw-related landforms and the topography they create are collectively termed as "thermokarst". However, the role of wildfires in developing thermokarst terrain remains quantitatively uncertain. Moreover, in comparison to the controlled warming experiments in Alaska , wildfires in arctic regions may also be viewed as uncontrolled disturbance experiments that aid in understanding the permafrost degradation processes.
Ice-rich permafrost deposits, known as the yedoma ice complex (yedoma), are widely distributed in the lowland of Alaska and Eastern Siberia . The greatest subsidence within the 2007 Anaktuvuk River tundra fire scar was identified in the yedoma upland by LiDAR . Yedoma is a unique permafrost deposit in terms of its extraordinarily high volume of ice (50-90 %) and organic-rich sediments. While the organic carbon trapped in permafrost regions is estimated to be twice that in the current atmosphere, permafrost thawing and related thermokarst processes may release the carbon as greenhouse gasses (CO 2 and CH 4 ) via microbial breakdown, which may further promote global warming . Thus, in order to estimate the volume of greenhouse gasses released, it is important to evaluate the volume of thawed ice associated with thermokarst processes in yedoma-rich areas.
Near the village of Batagay, Sakha Republic, Eastern Siberia (Figure 1), there exists the Batagaika megaslump, known as the world's largest retrogressive thaw slump, exposing roughly 50-90 m thick yedoma deposits on the north-east facing slope (e.g., . Thaw slumps are characterized by a steep headwall surrounding a slump floor and develop as a result of rapid permafrost thawing. The Batagaika megaslump was initiated at the end of the 1970s by deforestation but still appears to be growing . Considering this feature, it is worth considering whether new disturbances in the proximity will result in the formation of similar landforms. A wildfire incident occurred in July 2014 near Batagay, which, like deforestation, will change the ground thermal regime. Therefore, it is important to examine whether future catastrophic thermokarst development could be similarly initiated at the fire scar, whose area is much larger than the Batagaika megaslump (Fig 1b).
The first objective of this study was to assess the effectiveness of satellite Interferometric Synthetic Aperture Radar (InSAR) in detecting surface deformation signals due to wildfireinduced thermokarst over different temporal scales. InSAR has been used to detect long-term and seasonal displacements over several thaw-related landforms in permafrost areas (e.g., Iwahana et al., 2016;Antonova et al., 2018;. Although subsidence signals as a result of thermokarst associated with Alaskan wildfires have been detected using InSAR Iwahana et al., , 2016b, no such studies have been conducted on Siberian fires, to our knowledge. Also, all previous InSAR-based post-wildfire deformation mapping has been performed over relatively flat terrains, but no reports over hillslopes have been shown. Moreover, in contrast to previous studies, we employed two independent SAR imageries with distinct carrier frequencies and polarizations, L-band (1.2 GHz) HH-and C-band (5.4 GHz) VV-polarized microwave. Because the imaging geometries were different and had different sensitivities to the 3D displacement vector, we could not only take advantage of the performance of each sensor in mapping deformation signals but could also cross-validate the measurements by two InSAR data sets. 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 . 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  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 , 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 . However, a physical understanding of frost heave mechanisms has been established only during recent decades (e.g., . 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., . 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., . 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., . Although there is still an ongoing debate on the theory , we applied the simple, physics-based 1D theory of  to the observed frost heave signal so that we could physically interpret and explain the observed signals using reasonable parameters.

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 km 2 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 . 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  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 . 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 . 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 km 2 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 , 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 . 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.

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 lineof-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., Zwieback et al., 2015Zwieback et al., , 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 splitspectrum method . 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" , 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.  While ALOS2 has only imaged the area since 2014 its data acquisition interval is much longer than that of Sentinel-1 (Figure 2). Previous studies demonstrate that it is not possible to infer the total subsidence using pre-and post-wildfire SAR images, as the drastic changes in land cover cause low interferometric coherence . Additionally, JAXA changed the carrier frequency of PALSAR-2 in June 2015 ( Figure 2). Hence, monitoring long-term deformation using ALOS2 InSAR is possible only since October 2015. Conversely, frequent data acquisition in Sentinel-1 started only in 2017. Thus, we first performed an inter-comparison between ALOS2 and Sentinel-1 InSAR, focusing on the seasonal changes in 2017. We stacked Sentinel-1 interferograms to set the temporal coverages to nearly identical with those of ALOS2 (Figure 2). Stacking was necessary because we failed to derive long-term Sentinel-1 interferograms, as the burned areas quickly lost coherence. During the temporal interval of ALOS2 images Sentinel-1 had more cycles. Therefore, the number of Sentinel-1 stacks varied from three to eight (Figure 2).
Although L-band SAR is known to have better interferometric coherence than C-band SAR (e.g.,  our study indicated that Sentinel-1 could maintain a comparable interferometric coherence with L-band ALOS2 even during winter season. This is likely due to the short acquisition period of 12 days as well as the somewhat drier snow in the area in winter that allows microwave to reach the ground. A dry snow cover of depth less than 1 m is undetectable to microwave radiation, whereas over wet snow surface scattering dominates . The frequent data acquisition of Sentinel-1 since 2017 allowed us to examine detailed seasonal changes in surface deformation (Figure 2). Some Sentinel-1 InSAR pairs in early summer, however, did not show good coherence, possibly due to snow wetness.
In order to infer long-term temporal changes and cumulative displacements, we performed SBAS (Small Baseline Subset)-type time-series analysis Schmidt and Bürgmann, 2003), using 50 high-quality ALOS2 interferograms that included one-year-as well as short-term interferograms (Figure 3). We could estimate the average LOS-change rates between each acquisition epoch without assuming any temporal change models. In contrast to the original SBAS approach, we did not estimate DEM errors because the wellcontrolled orbit, as well as the precise TanDEM-X DEM, have no sensitivities to those errors.
In order to estimate the errors of the derived time series, we assumed each original SAR scene contained 0.2 cm errors, and made InSAR data covariance matrix, following the method of . The errors are relatively smaller than those in previous studies of SBAS analysis (e.g. 0.4 cm in 0.75 cm in Biggs et al., 2007), because, as noted earlier, we took out the long-wavelength phase trend form each InSAR image, and the analysis area is smaller (12 ×12 km) than previous studies.

Multispectral remote sensing of burn severity
Normalized burn ratio (NBR) is a useful multispectral remote sensing index to assess the impact of wildfire on vegetation. Vegetation reflects more strongly in the near-infrared (NIR) than in the shortwave infrared (SWIR) region, while a fire scar reflects more strongly in the SWIR. Utilizing this property, NBR is defined as NBR=(NIR−SWIR)/(NIR+SWIR). The difference NBR (dNBR), i.e., the difference between prefire NBR and postfire NBR, indicates burn severity (Key and Benson, 2006;Miller and Thode, 2007). Generally, when dNBR is greater than 0.66 the fire is regarded as "highly severe". We computed dNBR for the 2014 fire using Landsat 8, Band 5 (850-880 nm) and Band 7 (2110-2290 nm) images for near-infrared and shortwave-infrared, respectively, to associate the inferred subsidence distribution with burn severity.

One dimensional frost-heave theory based on premelting dynamics
We used the one-dimensional frost-heave theory as a tool to interpret the observed uplift signals. Inspired by one-way frost heave experiments ,  and  derived a steadystate heave rate V l of an ice lens, considering the force balance among thermo-molecular force F T , hydrodynamic force F μ , and overburden force F O (pressure P 0 ). Here, we assumed a constant heave rate V l , which may not necessarily reflect the actual observations shown below as well as in . However, this assumption simplified the theory, and we assumed that the observed heave rate did not change drastically over time.  proposed a non-dimensional heave rate v l of an ice lens as a function of its boundary position ξ l given: where μ, k 0 , and ρ are the viscosity of water, the permeability of ice-free soil, density of water, respectively. The quantity G ≡ ( L/ T m ) ⟨∇ T ⟩ has the same dimension as gravity and indicates thermo-molecular force when multiplied by the mass of displaced ice; L is the latent heat of fusion and T m is the bulk melting temperature. The first and second term in the bracketed numerator are proportional to F T and F O , respectively, while the bracketed denominator is proportional to F μ . The integral is performed along ξ ≡ z / z f , where z f is the position above (below) where ice saturation S s becomes non-zero (zero); z h indicates the position where hydrostatic pressure is achieved, and ϕ is the porosity of soil. The normalized overburden pressure and permeability are defined as p 0 ≡ P 0 / ρG z f and k ≡ k / k 0 ≥ 1, respectively.

Results
We performed an inter-comparison of ALOS2/Sentinel-1 interferograms, focusing on the seasonal changes in surface deformation. We then showed short-term deformation derived by Sentinel-1 and long-term deformation derived by time-series analysis of ALOS-2. Subsequently, we estimated the total volume of thawed excess ice. Although both satellite images covered the Batagaika megaslump we did not observe clear LOS changes as detected at the fire scar, which could be due to the lack of spatial resolution of the InSAR images. We compared the ALOS2 and stacked Sentinel-1 interferograms for five periods (Figure 4) and assessed their differences ( Figure 5). Despite differences in look directions both ALOS2 and Sentinel-1 indicated extensions in the LOS during periods (a, f) from the middle of June to the end of July and (b, g) from the end of July to the early October. Also, their deformation areas and amplitude were mostly consistent, suggesting that LOS changes were largely due to summer subsidence (see section 5.1 below for details). In terms of the spatial distribution of deformation signals, we noticed that the LOS changes over higher-elevation areas such as ridge and peak were insignificant, whereas the boundaries between the burned and unburned areas were clear. The north-western area, however, showed few LOS changes (see section 5.2 for the relationship between LOS changes and burn severity). During the period (c, h) from early October to early December both ALOS2 and Sentinel-1 indicated shortening in the LOS by an approximate 5 cm maximum, and the deformation areas and amplitude were quite similar. This observation presumably indicated frost heave in the early freezing period. In view of the previous three periods, both subsiding and uplifting areas were nearly the same.

Seasonal deformation and comparison of ALOS2/Sentinel-1 interferograms
The following period (d, i) from early December to the middle of March also included the winter season with much colder air temperatures, but we did not observe any significant deformation signals, indicating that frost-heave virtually stopped in early December. While the good interferometric coherence during mid-winter was an unexpected result, we speculate that it could have been due to drier, lower amounts of snowfall.
In the periods (e, j) from the middle of March to early June, both ALOS2 and Sentinel-1 suffered from decorrelation, and we could not identify clear deformation signals. However, in light of Figure 6 below, each of the Sentinel-1 interferograms had overall good coherence with the exception of the data acquired in the middle of May. These observations suggested that the decorrelation may be attributable to the rapid changes on the ground surface during the initiation of the thawing season when the air temperature rises above the freezing point and the active layer begins to thaw.  (Figure 4e and 4j), we could not estimate differences due to coherence loss. Figure 5 shows the differences between ALOS2 and Sentinel-1 InSAR data with nearly identical periods, which may help in cross-validating the measurements and understanding the actual deformation processes. The estimated differences and their 2σ scatter were 0.5±1.2cm (Fig 5a),0.7±2.3cm (Fig 5b),0.3±1.3cm (Fig 5c), and 0.6±0.3cm (Fig 5d), with mean of 0.5±1.5 cm. The differences and their variances were variable over time but apparently indicated some systematic trends. For instance, over the east-facing slopes, the differences were almost always positive (This is discussed more comprehensively in section 5.1).
The Sentinel-1 interferograms for 2017 demonstrate that the progress of deformation was not at a constant rate ( Figure 6). The most rapid deformation took place in June (periods 1 and 2) with no substantial deformation in July (period 3) and started to subside again in August (periods 4-6). We found that the subsidence occurred sporadically over time and space and that the burned area did not uniformly subside. For periods 4, 5 and 9 we were unable to perform phase unwrapping at specific locations near the ridge and the boundaries between the burned and unburned areas. These unwrapping errors were responsible for the localized, large differences observed in Figure 5b. We confirmed the presence of low coherence bands along the unwrapping errors, which may suggest large phase jumps due to large displacements during the 12 days; enigmatically, no such line-shaped low coherence was detected in the long-term ALOS2 interferograms. Moreover, Figure 6 demonstrates that the frost heave started in late September, which was missed in the periods (b) and (g) of Figure 4, and that the absence of any deformation signals lasted from early December to May of the following year. We will physically interpret the absence of deformation signals during the coldest season in section 5.3.  Table 1; imaging was performed by ascending, right-looking orbit. Warm and cold colors indicate LOS changes away from and toward the satellite, respectively. Black dashed line indicates the boundary between the burned and unburned area confirmed with Landsat optical images.
Figures 7a-7h show ALOS2 interferograms, each of which covers nearly one-year after October 2015 with some overlaps in its temporal coverages. Figure 7a, derived at the earliest period after the fire, indicates the maximum one-year subsidence to be as much as 10 cm or more.
If the amplitude and timing of seasonal subsidence/uplift cycle are invariable over time, a one-year interferogram will tell us only the irreversible displacements regardless of the acquisition times of master/slave images, which corresponds to the "pure ice" model in . Figure 7 sequentially shows the periods from October 2015 to June 2019 and indicates that the yearly subsidence rate slowed down. However, the variations of the oneyear LOS changes in Figures 7 suggest that the actual deformation processes were more complex. Figure 8a shows the cumulative LOS changes from October 2015 to June 2019 derived from SBAS-type time-series analysis, and that the maximum LOS extension reached as much as 25 cm; the 2σ errors for Figure 8a were ±1.5cm. Considering that the LOS changes during the first year after the 2014 fire were not included, the total LOS changes were presumably much greater than 25 cm, which meant that the subsidence was greater than 30 cm on account of the 36° incidence angle. As mentioned earlier, however, the higher-elevation areas such as the ridge did not undergo significant deformation, which probably would have been the case even during the first year after the fire. In addition to the high elevation areas, we realized clear contrasts in the LOS changes between the east-and the west-facing slopes near the northwestern area and the central north-south trending ridge; this spatial heterogeneity could also be recognized in Sentinel-1 (Figure 4). Their possible mechanisms comparing the burn severity ( Figure 8b) and local landform (Figure 8c) are discussed in section 5.2.
We show the estimated time-series data at four representative sites (Figures 9a-9d), whose locations are indicated in Figure 8a. The sites (a) and (b) underwent nearly the same cumulative LOS changes by roughly 20 cm but were located at different slopes that are 4.3 km apart. On the other hand, the cumulative LOS changes at the site (d) were relatively small (approximately 10 cm). The site (c) located in the ridge did not show either significant seasonal or long-term deformation.

Estimating the total volume of thawed excess ice
Post-wildfire deformation over a permafrost area presumably consists of two contributions: (1) irreversible subsidence due to melting of ice-rich permafrost below the active layer, and (2) seasonally cyclic subsidence and uplift due to freeze-thaw of the active layer . In order to separate the two processes from the observed deformation data,  used independent ground-measured ALT data to predict the ALT contribution to total subsidence. Ground-measured pre-fire ALT data were not available at this study site. Given the temporal evolution of post-wildfire deformation data (Figures 9a-9d), however, we may regard the cumulative deformation in Figure 8a as being due to irreversible subsidence during the period between October 2015 and June 2019, and estimate the total thawed volume as 3.56 ± 2.24 ×10 6 m 3 ; the error bar is based on the root mean square of the no-deformation signals outside the burned area, which is multiplied by the burned area. However, in view of the temporal evolution in Figure 9, we could speculate that a much larger deformation was also taking place immediately after the 2014 fire until October 2015, during which, unfortunately, no deformation data are available. Thus, this estimate should be viewed as a lower estimate, with the actual volume of thawed permafrost possibly being much greater.
Nevertheless, despite its much smaller area size (Figure 1b), the thawed volume at the Batagaika megaslump is greater than 2.5 ×10 7 m 3 , an order-ofmagnitude larger than our estimate above. Moreover, the thaw-subsidence rate at the fire scar is slowing down (Figure 9). We discuss the possibility of the another megaslump emergence at the fire scar in section 5.2.

Similarities and differences in the ALOS-2 and Sentinel-1 interferograms: implication for insignificant slope-parallel sliding
Taking into account the imaging geometries of the ALOS2 and Sentinel-1, we could comprehensively interpret the differences in Figure 5 and also infer the actual deformation processes. The weights multiplied to the 3D-displacements, (U ew , U ns , U ud ), to compute LOS changes were +0.573, +0.132, and −0.809 for ALOS2 and −0.583, +0.236, and −0.777 for Sentinel-1, respectively; eastward, northward and upward displacements were taken to be positive. Assuming the LOS changes of the two sensors are identical (which is roughly the case in Figure 4), and no north-south displacement U ns , the constraint on the east-west and up-down displacements can be derived as U ew :U ud = 0.032:1.156. The assumption of zero U ns might appear unrealistic but can be reasonable over the east-and west-facing slopes, which incidentally cover a broad area of the fire scar. As this constraint indicates the dominance of vertical displacement, we can infer that slope-parallel sliding did not take place over the eastand west-facing slopes.
In the thawing season when the vertical displacement is downward (negative), the previous constraint on the two displacements also indicates that the east-west displacement should always be westward (negative), regardless of the slope. As this is physically implausible, we may assume that east-west displacements were virtually zero over both the east-and westfacing slopes. We can thus infer a pure vertical subsidence without any east-west displacements during the thaw season. Hence, the differences between ALOS2 and Sentinel-1 in the thawing seasons (Figures 5a and 5b) will be simply equal to −0.032U ud . Therefore, we can expect systematically positive differences in the thawing season, regardless of the eastand west-facing slopes, which appear consistent with observations (Figures 5a and 5b). Quantitatively, however, the mean differences of 0.5-0.7 cm are too large to be attributable to the geometric difference alone, on account of the subsidence by as much as 5 cm or more. Here, we hypothesize the possible impact of soil-moisture changes, which can reach ~10 % of the carrier wavelength (Zwieback et al., 2015(Zwieback et al., , 2016. As changes in soil moisture generate larger apparent LOS changes in L-band than in C-band InSAR, the observed differences can be likely. In contrast to thaw subsidence frost-heave is more likely to occur towards the slope normal direction. Assuming the magnitude of slope-normal uplift, U f , over a slope with gradient θ, the differences between ALOS2 and Sentinel-1 would be U f (1.156 sinθ −0.032cosθ) assuming zero U ns . We estimated |θ|=1.58°, which corresponds to 55 m height difference over 2 km horizontal distance and was fairly consistent with the slope of the studied area. Meanwhile, the differences can also be considered 1.156U ew −0.032U ud , which indicates additional positive and negative effects on the east-and west-facing slope, respectively. Indeed, Figure 5c appears to depict clearer contrasts in sign on the east-and west-facing slopes. Moreover, the impact of changes in soil moisture are likely much smaller in the colder season than in the thaw season, which may explain the smaller differences in freezing seasons.
As we derived the differences over the no-deformation season (Figure 5d), we can attribute them to the atmospheric effect on ALOS2 and Sentinel-1 interferogram (Figures 4d and 4i). The overall positive differences are likely because the spatial scale of atmospheric delay was greater than the fire scar area.
Previous reports of thermokarst subsidence after fire have focused on relatively flat areas as those at the 2002 tundra fire in the central Seward Peninsula, Alaska (Iwahana et al., 2016b), the 2007 Anaktuvuk River tundra fire , and the 2009 Big Creek Fire in the Alaskan Yukon River basin . As such, in addition to the broad subsidence detected by InSAR, polygonal patterns associated with ice wedge degradation became clearly visible 4-7 years after the fire by highresolution optical and LiDAR remote sensing Iwahana et al., 2016b). At the studied hillslopes, in contrast, no such polygonal patterns are likely to be detected. Nonetheless, the dominance of vertical displacements with little slope-parallel sliding indicate that rapid active-layer detachment sliding (ALDS) events were insignificant. In contrast, many ALDS events triggered by fire have been mapped at Mackenzie Valley, Canada, whose length could sometimes reach hundreds of meters . If ALDS event with such length occurred, we could have observed significant loss of interferometric coherence. It is possible, however, that local ALDS events occurred but were undetected because of the coarse resolution (~10 m) of InSAR images. Because the subsidence was caused by thawing of ice-rich permafrost, meltwater should have been supplied at the base of active layer. Considering the mechanisms of ALDS , porewater pressure increase might have been not enough to reduce the effective overburden stress and to initiate significant slope-parallel sliding. This is possibly because the meltwater could have drained through the gullies. However, in view of the significant uplift signals over the burned area even years after the fire, the meltwater is still likely to be undrained on the slope. If there were further enough water input by, for instance, warmer days and/or heavy rain, significant ALDS events may take place in the future.

What controls the heterogeneous distribution of subsidence magnitude? Possible emergence of another megaslump
The cumulative subsidence magnitude was spatially variable but showed some systematic changes. In addition to the ridges and peaks the west-facing slopes showed significantly smaller subsidence than the east-facing slopes (Figures 4 and 8a). To interpret the spatially heterogeneous subsidence, we associated burn severity and local landform with the cumulative subsidence ( Figure 8). In light of the inferred dNBR (Figure 8b), which ranged from 0.2 to 0.4, the burn severity was moderate rather than high. Also, the burn severities were spatially less heterogeneous than those of cumulative subsidence and local landform. In fact, we could even identify deformation-free areas having even higher burn severity. Thus, although the fire undoubtedly initiated the subsidence, the burn severity did not control the subsequent cumulative magnitude .   37   519  520  521  522  523  524  525  526  527  528  529  530  531  532  533  534  535  536  537  538  539  540  541   542  543   544  545  546  547  548  549  550  551  552  553 Notably, however, gullies were clearly more developed on the east-facing slopes than on the west-facing slopes (Figure 8c), which were confirmed to be present at least back in 1991 by Landsat image. Considering the striking correlation between the development of gullies and the larger subsidence, there is high likelihood of a causal relationship between them. Similar dependence on the slope aspect was reported by Lacelle et al (2010,2015), who found that hillslope thaw slumps in the Richardson Mountains-Peel Plateau, northwest Canada, predominantly developed on the east-facing slope.  interpreted that the active layer on the east-facing slope might be thinner because of lower amount of insolation than on the south-and west-facing slopes, which promoted a triggering mechanisms of thaw slumps because the ice-rich permafrost was closer to the surface. Although the broadly subsiding areas are not so-called thaw slumps, thinner active layers on the east-facing slope are likely and can consistently explain both the larger subsidence and the rich development of gullies. This hypothesis can be tested either by examining the surface deformations at the 2018 and 2019 fire scars and other fire scars across Siberia and other boreal regions or by performing field-based thaw-depths measurement.
The recent slowdown of the subsidence rate ( Figure 9) may suggest that the 2014 fire scar could stabilize in the near future. However, although it depends on how quickly the vegetation is recovered, we do not preclude the possible emergence of another megaslump particularly on the east-facing slopes. In order to initiate thaw slumps, ice-rich permafrost needs to be exposed at the surface , at which the initial headwall and slump floor are formed. In contrast to the thaw slumps near shorelines, coastlines and riverbanks (e.g., Kokelj et al., 2009), no mechanical erosions by waves and currents are effective on hillslopes like the studied area. For the development of retrogressive thaw slumps (RTS) on hillslopes, Lacelle et al. (2010) suggested that ALDS triggered by meteorological events could remove the overlying active layer and expose the ice-rich permafrost. Although no large-scale ALDS events were detected during the studied period, they might take place as discussed in the previous section. Moreover, Figure 8a indicates that the subsidence magnitude becomes larger toward upslope, and there are clear boundaries between the subsiding and non-subsiding portions, where an initial headwall for RTS could be exposed. Once an initial headwall has formed, subsequent retreat rate is rapid on the order of several meters per year . Thus, in order to monitor the early formation process of RTS in more detail, long-term radar remote sensing with higher spatial and temporal resolution would be necessary and promising. 39   554  555  556  557  558  559  560  561  562  563  564  565  566  567  568   569  570  571  572  573  574  575  576  577  578  579  580  581  582  583  584  585 586 587 Figure 10. Non-dimensional heave rate profiles of an ice lens as a function of its boundary position, based on the analytical model by . Five cases of nondimensional overburden pressure p 0 and porosity ϕ are shown.

Interpretation of frost-heave signals
In order to interpret the detected frost heave signal on the basis of the theory by , we first examine the sensitivity of the heave rate on the normalized overburden pressure p 0 and porosity ϕ. Figure 10 shows five cases of non-dimensional heave rate profiles as a function of the ice lens boundary position ξ l , indicating that the maximum heave rate is mainly controlled by the normalized overburden pressure p 0 and is somewhat insensitive to the porosity ϕ. Details of the heave rate profiles will depend on the assumed models of permeability and ice saturation, but the qualitative characteristics are not altered . There exist two positions that give the same heave rate, but only the branch with smaller ξ l is stable .
We can attribute the clear contrast in the frost heave signals inside and outside the burned area to the differences in the normalized overburden pressure p 0 . Because the mechanical overburden pressure P 0 will not significantly differ from the inside to the outside of the burned area, the larger frost heave rate in the burned area would be caused by larger temperature gradient G and/or deeper frozen depth z f . Owing to the removal of vegetations and surface organic layers over the burned area, the larger temperature gradient G than that of the unburned area is likely more marked in the early freezing season and may generate a greater thermomolecular force that will effectively reduce the normalized overburden pressure. We may also interpret the absence of frost heave signals in mid-winter as due, probably, to the smaller temperature gradient G than that in late fall/early winter; if frost heave were controlled by temperature instead of temperature gradient, we would expect even more significant signals during the much colder part of the season. The deeper frozen depth z f is also likely due to the loss of surface vegetation and should supply more water for frost heave. From the end of September to the middle of November 2017, Figure 6 shows LOS changes by approximately 1.5 cm over 12 days toward the satellite that corresponds to an approximate 1.9 cm uplift. Assuming a constant-rate frost heave, this corresponds to a heave rate of 1.8 × 10 -8 (m/s). The most critical parameter controlling heave rate is the permeability for ice-free soil k 0 , which can vary by orders-of-magnitude, while other parameters are well-constrained. We may fit our observed heave rate with the ice-free permeability, k 0~1 0 -17 (m 2 ), which is a likely value in view of the three cases in .
Here we comment on the modeling of uplift signals as caused by in situ freezing of pore water into ice . The in situ freezing model is simple, and can explain the timing, duration, and magnitude of uplift signals, if one assumes such pore water in the active layer. However, because the Stefan function approach in  is essentially controlled by atmospheric (or ground) temperature that is rather homogeneous over this spatial scale, it is difficult to account for the observed heterogeneous distribution of uplift signals. The distribution of uplift signals was closely correlated with that of subsiding signals, which led us to interpret that the permafrost thaw and its incomplete drainage could become a water reservoir for ice lens formation and frost-heave. The frozen pore ice within the soil and the ice lens formed by water migration are totally different in terms of their formation mechanisms and subsequent forms of ice. From a geomorphological perspective, the presence of ice lenses will play a role in reducing the strength of soil and potentially initiating ALDS, because porewater pressure will increase at the front of thawing, whereas pore ice within the soil would simply stay as pore water with little impact on the landform.
We also recognize, however, that the microphysics-based theory adopted in this study is developed in 1-D geometry and is based on the assumption of "frozen fringe", a region where liquid freezes into ice through the pores of soil. Some laboratory experiments did not support the presence of frozen fringe (e.g., , and the "fringe free" frost heave theory has also been proposed; see  for review. In addition to the controlled lab experiments and theoretical developments, more detailed observations of natural frost heave signals are becoming possible and might help better understand the physics of frost heave and its geomorphological consequences.

Conclusions
We used L-band and C-band InSAR to detect post-wildfire ground deformation at Batagay in Sakha Republic, showing not only subsidence signal during the thawing season, but also uplift during the early freezing season and virtually no deformation in midwinter without loss of coherence. Time series analysis allowed us to estimate cumulative displacements and their temporal evolution, as quality interferograms could be obtained even in the winter season. We found that the thawing of permafrost in the burned area lasted three years after the fire, but apparently slowed down after five years. During the studied period, no significant slopeparallel sliding was detected, and the post-wildfire deformation was mostly subsidence. Despite the rather homogeneous burn severity, the cumulative subsidence magnitude was larger on the east-facing slopes and showed a clear correlation with the development of gullies, suggesting that the east-facing active layers might have been originally thinner. Shortterm interferograms (2017Shortterm interferograms ( -2018 indicated that the subsidence and uplift was clearly enhanced compared with the unburned site. We have thus interpreted the frost heave signals within a framework of premelting dynamics. Post-wildfire areas are a focus of permafrost degradation in the Arctic region.