Surface displacement revealed by L-band InSAR analysis in the Mayya area, Central Yakutia, underlain by continuous permafrost

Recent increases in global temperature have stimulated permafrost degradation associated with landform deformation caused by the melting of excess ground ice (thermokarst). Central Yakutia is underlain by ice-rich continuous permafrost, and there are complicated permafrost-related features in forested and deforested areas. This situation makes thermokarst monitoring necessary over a wide area to achieve a better understanding of its dynamics. As a case study, we applied L-band InSAR analysis to map surface subsidence due to thermokarst in this area and to demonstrate the suitability of L-band SAR for such monitoring. Our results show that InSAR detected subsidence/uplift signals in deforested areas and alasses; whereas, there were few ground deformation signals in forested areas with middle coherence. The InSAR stacking process, including both seasonal and inter-annual displacements, showed subsidence in deforested areas during 2007–2010 and 2015–2018, in the range of 0.5–3 cm yr−1. We also estimated the inter-annual subsidence to be up to 2 cm yr−1 during 2015–2018, using InSAR pairs that spanned the same seasonal interval but in different years. The magnitude of subsidence and the spatial patterns are qualitatively reasonable as thermokarst subsidence compared to observations using field surveys and high-resolution optical images. L-band InSAR was effective in maintaining coherence over a long period for a partially forested thermokarst-affected area, which resulted in deriving the inter-annual subsidence by the stacking using four interferograms. The advantage of the persistent coherence in L-band InSAR is crucial to better understand thermokarst processes in permafrost regions.


Introduction
Thermokarst is a characteristic landform that results from the thawing of ice-rich permafrost and melting of massive underground ice (van Everdingen 2005;French 2017). Thermokarst processes can form large depressions in the ground, leading to ground inundation. This irreversible subsidence leads to the formation of thermokarst lakes and changes the entire landscape. Thermokarst-induced geomorphological and hydro-ecological changes have been reported in several studies through in situ observations, satellite images, and model simulations (e.g., Zakharova et al. 2015;Ulrich et al. 2017;Aas et al. 2019). Thermokarst has also been found to damage or destroy infrastructure (Shiklomanov et al. 2017;Hjort et al. 2018). To mediate the damage and adapt to landform changes, understanding about the deformation rates and spatial distribution of the phenomenon on a regional scale is essential. It is necessary to obtain surface displacement data for a better understanding of the thermokarst dynamics.
In Eastern Siberia, specifically in Central Yakutia ( Fig. 1), ground ice constitutes 50-80% by volume (Brouchkov et al. 2004). The mean annual temperature in Yakutsk is − 9.7 °C, which increased in Central Yakutia Open Access *Correspondence: abe.takahiro2@jaxa.jp 1 Earth Observation Research Center, Japan Aerospace Exploration Agency, Tsukuba 305-8505, Japan Full list of author information is available at the end of the article in the early 1990s and 2005(Fedorov et al. 2014, during which thermokarst subsidence was activated Konstantinov 2003, 2008). Short-term dynamics have been examined in several studies (e.g., Brouchkov et al. 2004;Fedorov and Konstantinov 2009). In Yukechi, close to Mayya (Fig. 1), surface dynamics have been monitored by a research group at the Melnikov Permafrost Institute since 1992, where the annual subsidence of up to 8 cm was measured on a flat area of interalas meadows (Fedorov and Konstantinov 2009); an alas is the final geomorphological stage of old thermokarst development (e.g., van Everdingen 2005). Young thermokarst lakes in Yukechi increased fourfold from 1980 to 2012, and new lakes have been forming (Fedorov et al. 2014). The differences in the behavior and distribution of the lakes in Yukechi have been studied using satellite optical images by Ulrich et al. (2017). Saito et al. (2018) recently examined thermokarst landforms using an Unmanned Aerial Vehicle and the Structure from Motion technique on a disused airfield and on arable land in the Churapcha area. They revealed detailed distributions of high-centered polygons with average diameters of 11.6 and 7.4 m, respectively, and average subsidence rates of 2.1 and 3.9 cm yr −1 since 1990 were estimated.
Central Yakutia is located in a continuous permafrost region (Brown et al. 2002) that exhibits complicated permafrost-related features. Many alasses developed thousands of years ago, and other open areas experienced surface disturbances by deforestation for land use in the 1970s. Under these circumstances, thermokarst development has been ongoing in both forested and deforested areas for several decades. Field observations were conducted to determine the rate of thermokarst-related subsidence by point measurements as mentioned above. Satellite and airborne optical images have enabled better understanding of the development of thermokarst lakes, but they cannot derive the rate of surface displacement due to thermokarst settlement. Therefore, both techniques are limited in their application to monitoring thermokarst in the entire Central Yakutia.
Synthetic aperture radar (SAR) interferometry (InSAR) is a processing technique for detecting surface displacements that has been widely used in applications such as crustal deformation (e.g., Massonnet et al. 1993) and glacier motion (e.g., Joughin et al. 2010). SAR data have been recently used to monitor permafrost's related deformations in Alaska (Liu et al. 2010(Liu et al. , 2015Iwahana et al. 2016a), Canada (Short et al. 2011(Short et al. , 2014, Lena Delta Chen et al. 2018a;Strozzi et al. 2018), and Tibet (Chen et al. 2018b). Liu et al. (2015) used the Advanced Land Observing Satellite/Phased Arraytype L-band Synthetic Aperture Radar (ALOS/PALSAR) InSAR to reveal the thermokarst settlement near Deadhorse along the North Slope, Alaska. They demonstrated that InSAR is an effective tool for mapping and studying the thermokarst process. However, InSAR has not been used to assess deformation related to permafrost in Central Yakutia; therefore, we need to examine the usefulness of InSAR in this area.
Synthetic aperture radar satellites carry SAR sensors using different microwave wavelengths. TerraSAR-X carries an X-band SAR (wavelength 3.1 cm) that generates high-resolution (1-3 m) images over a short interval (11 days). The Sentinel-1A/B constellation has enabled us to generate a large number of interferograms with a revisit interval of 6 or 12 days, and is available to end users at no cost. The available C-band data of Sentinel-1 have high potential to allow understanding seasonal thaw dynamics. Recent papers on permafrost studies using Sentinel-1 data emphasize its usefulness in deriving seasonal displacement (e.g., Chen et al. 2018a;Strozzi et al. 2018). Although the short microwave wavelengths used by the X-and C-band SAR sensors are more sensitive to smaller surface displacements than longer wavelengths such as L-band SAR, the disadvantage is that the coherence in X-and C-band data quickly decreases with time separation; this makes it difficult to generate useful, coherent interferograms for intervals of more than a month. Therefore, examining the inter-annual displacement related to the thermokarst is particularly challenging. Wang et al. (2017) evaluated the influences and conditions for applying InSAR using a number of X-and L-band interferograms in a discontinuous permafrost region of northern Quebec, Canada. They reported good L-band results with deeper penetration through vegetation that maintained coherence over time and space, indicating that L-band SAR has a notable advantage over X-band SAR in monitoring displacement. Many studies have targeted ground-surface displacement in periglacial landscapes using ALOS/PALSAR L-band InSAR (e.g., Short et al. 2011;Chen et al. 2013Chen et al. , 2018bJia et al. 2017;Dini et al. 2019;Michaelides et al. 2019). However, there have been a limited number of studies on the inter-annual measurements of thermokarst subsidence by L-band InSAR. ALOS and ALOS-2 operated by Japan Aerospace Exploration Agency (JAXA) have been the only L-band SAR satellites in the last two decades, and L-band InSAR has great potential for revealing surface displacement to better understand the thermokarst process.
The aim of this study is to apply L-band InSAR for Central Yakutia and to measure inter-annual subsidence due to thermokarst, supported by ground surveys. We focused on the spatial and temporal changes in surface subsidence due to thermokarst in deforested areas cultivated a few decades ago, and in alasses developed thousands of years ago, and interpreted the signals detected. We also assessed the utility and limitations of L-band InSAR as a monitoring tool for thermokarst through observation of the measurements. This paper presents a case study of L-band InSAR utilization for thermokarst monitoring.

Study area
Mayya (Fig. 2a) is a rural locality located on the right bank of the Lena River, and 40 km southeast of Yakutsk. Mayya is mostly surrounded by forest (Fig. 2b), although the town was deforested for farming, primarily in the 1970s (Fig. 2a, b). Mayya is located on fairly flat ground with an elevation of ~ 140 m, and the southern part of the town is higher with an elevation up to 200 m (Fig. 2c). Alasses can be identified by their slightly lower elevation and greenness compared to the surrounding yellow interalas meadow (Fig. 2c). The long-term average thickness of the active layer of thermokarst-affected areas ranges from 2 to 3.5 m, while those in intact areas range from 1.8 to 2.7 m in the Mayya region. Mayya is representative of residential areas where thermokarst development has been reported in Central Yakutia, but there are no studies on permafrost's related surface deformation. It is essential to know where and how much the surface subsidence occurs due to thermokarst for a better understanding of the current situation in Mayya.

SAR data processing
In this study, we used L-band SAR data obtained by ALOS/PALSAR and the subsequent satellite sensor ALOS-2/PALSAR-2. ALOS was launched in January 2006 and terminated in May 2011. ALOS-2 was launched in May 2014 and is still in operation as of February 2020. We selected datasets from the stripmap mode in ALOS and ALOS-2 because of the higher observation frequency and spatial resolution than those of the other modes.
From 2007 to 2010, ALOS/PALSAR acquired 17 images in the stripmap mode in both fine beam single polarization (FBS) and fine beam double polarization (FBD) over Mayya (path 409, frame 1230). FBS only obtains HH polarization images, while FBD obtains a dual polarization image of HH and HV. InSAR results using snow-covered SAR images showed strong decorrelation. The pairs with larger perpendicular baseline (B-perp; distance perpendicular to line-of-sight between the positions of the satellite at different times) contained unreliable deformation signals because large B-perp causes geometric decorrelation and topographic phase components in the presence of digital elevation model (DEM) errors. Thus, the ALOS data obtained between June and the beginning of October, and interferometric pairs with sufficiently high coherence (a threshold of 0.4) were chosen. As a consequence, we used four interferograms using four images from the ALOS/PALSAR FBD mode with 20-m resolution (Table 1).
From 2014 to 2018, ALOS-2/PALSAR-2 obtained seven scenes of the stripmap mode 3 (SM3) with 10-m resolution (hereafter SM3; path 126, frame 1230, beam F2_6) and six scenes of the stripmap mode 1 (SM1) with 3-m resolution (hereafter SM1; path 28, frame 2370, beam U2_7) over Mayya. The center frequency modification of the beam F2_6 (off-nadir angle is approximately 32°) in SM3 was carried out in June 2015 (JAXA 2015), which means interferometry between the F2_6 data acquired before and after June 1 2015 is undesirable because of the significant loss in range bandwidth (Natsuaki et al. 2016). Thus, we did not use the SM3 data acquired before June 1, 2015. The orbit of ALOS-2 has been controlled within a 500-m radius tube around the reference orbit (e.g., Ohki et al. 2018). Therefore, we expect that B-perp has little effect on InSAR results on the pairs, assuming no large DEM errors. Since the SM3 data were mostly acquired at the end of September or the beginning of October of each year, we chose only those interferometric pairs with sufficiently high mean coherence conditions (threshold of 0.3). One SM1 image was acquired in August 2016 and the others were obtained in August-September 2018. To derive seasonal displacement, we chose the SM1 data obtained in August-September 2018. Finally, we selected 11 interferograms using seven images of SM3, and six interferograms using four scenes of SM1 for ALOS-2/ PALSAR-2 (Table 1).
We used GAMMA software (Wegmüller and Werner 1997) to process the SAR data and generated Single Look Complex (SLC) data from Lv1.0 data in ALOS and Lv1.1 data in ALOS-2. This processing was described by Strozzi et al. (2018). After precise co-registration between two SLC images, we generated interferograms with a multilooking factor of 3 and 2 pixels in slant-range and of 9 and 4 pixels in azimuth for ALOS and ALOS-2, respectively.
The topographically related phase was removed using ALOS World 3D (AW3D), a 5-m high-resolution digital surface model derived from ALOS. We also applied Goldstein-Werner's adaptive spectral filter with an exponent of 0.6 to smooth the signals (Goldstein and Werner  (Costantini 1998), was then processed with 2 × 2 pixels of multi-looking factor and threshold values of coherence as 0.4 and 0.3 for ALOS and ALOS-2 to mask out the unreliable phase. The spatial resolution after terrain-corrected geocoding projecting onto the UTM coordinate was set to 25 m. Finally, the line-of sight (LOS, distance in satellite and ground) change detected by InSAR was converted to a vertical displacement by dividing by the cosine of the incidence angle (e.g., Liu et al. 2010;Iwahana et al. 2016a). We consider this assumption to be suitable because thermokarst-induced displacement is almost dominant in vertical deformation and because the effect from topographic slope is small. Some interferograms were affected by orbital errors, and ionospheric and atmospheric disturbances. The effects of orbital errors and ionospheric disturbances are sometimes shown in a linear long-wavelength trend. Our study area is spatially limited to ~ 10 × 7 km 2 and the subsidence caused by thermokarst is a localized deformation. Thus, we were able to model and subtract the large-scale trend from the original interferograms by fitting a 2D linear function. We decided that the reduction was sufficient, and no ionospheric correction, such as the Split Spectrum Method (e.g., Gomba et al. 2016), was necessary. The effects of the laying of troposphere sometimes cover an interferogram with scales of several kilometers, resulting in apparent displacements of up to a few centimeters or more in an interferogram (e.g., Hanssen 2001). Moreover, the intensity of the effect is often proportional to topography, especially in mountainous areas. On the other hand, our target area is situated in a relatively flat fluvial terrace of the Lena River, and the regional climate is very dry (annual precipitation is about 250 mm). Thus, we considered the tropospheric effect to be sufficiently small to be neglected. There could also be a turbulent atmospheric effect in an interferogram; these atmospheric effects are sometimes estimated from the variogram and used to obtain a geophysical model parameter through multiple optimizations (e.g., Sudhaus and Jónnson 2009).
We had to set a spatial reference point (i.e., ground control point) that is not displaced in interferograms. However, it was difficult to find a stable point in our study area because the ground can move through freezing and thawing in a permafrost region. Thus, we set a reference point based on these ideas; (I) The InSAR phase at a point with low coherence is unreliable as a stable point; (II) a point in an alas is considered to be more stable as a final stage of thermokarst than that in other ground areas underlain by massive ground ice in this region. Therefore, we calculated averaged coherence maps using FBD, SM3, and SM1 (Additional file 1: Fig. S1) and selected a reference point in an alas ('x' in Figs. 3,4 and 5). This assumption may cause errors related to the selection of the reference point .
To extract the averaged surface displacement rate, we performed a stacking analysis on interferograms according to Eq. (1), where Pair is each interferogram converted to displacement rate and tdays represents the total number of days including the thawing season (defined as the period from June 1 to September 30) for the two acquisitions. We assumed that surface subsidence occurs in the thawing season and no subsidence occurs during winter freezing in our datasets. Performing such weighting reflects the actual period that contributes to the signals of surface subsidence rather than using the total number of days for the two images. N is the number of interferograms. The resolution of the stacking results was 25 m, which is the same as that for each interferogram.

Field observations
To measure surface displacement, we located survey nails at 1-m intervals along two 30-m transects, at four plots (A-D in Fig. 2d) within a ~ 1 km 2 area located in an abandoned agricultural field (currently used as meadow, Fig. 2e). The relative heights of the ground surface adjacent to the nails, which were attached to the ground (1) InSAR stacking cm yr −1 = Pair 1 × tdays 1 + Pair 2 × tdays 2 + · · · Pair N × tdays N N 1 tdays N , surface, were determined by an automatic optical level (Leica NA720) with 0.5-cm vertical accuracy with respect to the tops of nearby permafrost borehole pipes. Different nearby borehole pipes were used for the pairs of Plots A/B and Plots C/D. Plot E was located in the center of an area deforested and abandoned in the 1990s (Fig. 2a). We set seven permanent survey markers within a 30-m line and repeatedly measured the relative heights of the ground surface adjacent to the markers against a control point in 2015 and 2018. Assuming that the control point had been stable, we calculated an average displacement rate over the 3-year interval. However, the control point marker at Plot E was not fixed in permafrost because there were no reliable benchmarks for the leveling survey and the plot may have subsided, leading to underestimates in the absolute subsidence at other markers. In this area, local people witnessed ground-surface subsidence after clear-cutting. The initial gentle slope landscape changed to depressed, bumpy relief. An inter-annual upheaval due to permafrost aggradation was less probable because Plot E has been experiencing thermokarst subsidence because of the surface disturbance. Assuming the temporal stability of the control point during our observation period, we treated the measured subsidence rate as a minimum value in 2015-2018.
At Plots A-D, we performed a number of leveling measurements along the transects, which were performed in late September in 2017 and 2018. Ground Fig. 3 Vertical displacement derived from unwrapped interferograms using a ALOS FBD, b ALOS-2 SM3, and c ALOS-2 SM1 data. The positive and negative values correspond to uplift and subsidence, respectively. White shows areas of no data due to low coherence and/or no data in AW3D. The 'x's indicate the reference points of the interferograms. The two areas marked A1 and A2 indicate the locations of the obvious signals (the legends are also used hereafter). The red arrows in (b) also show the obvious signals excluding those in (a) surface displacement was determined by the difference between the two measurements at a 1-year interval. The two times standard error of the difference was considered as the measurement error for each plot. The FBD interferograms obtained on July 23 and September 7, 2009 show two remarkable areas with subsidence signals (A1 and A2 in Fig. 3a). The magnitude is up to 3.5 cm for 46 days. Figure 3b shows the SM3 interferogram obtained on September 29, 2017 and September 28, 2018. The two prominent subsidence areas (A1 and A2) instead show a 2-3-cm uplift in Fig. 3b. On the other hand, marked subsidence signals of 1-2 cm (red arrows in Fig. 3b) were also measured. The SM1 interferogram (Fig. 3c) shows a more spatially detailed and coherent image in surface displacement than the previous two interferograms because of a higher spatial resolution (3 m) and a shorter revisit time (42 days). The signals in A1 and A2 are represented as subsidence of 2-3.5 cm, and other areas with subsidence signals can also be seen in Fig. 3c. Figure 4a shows the stacked vertical displacement rate derived from the FBD data for 2007-2010 (Table 1 and Fig. 4d). Compared with the result in Fig. 3a, the signals  Table 2. The red arrows in (a) and the areas surrounded by solid lines are obvious subsidence signals. The signals indicated from T1 to T5 in (c) show obvious inter-annual subsidence. d Temporal distribution of InSAR pairs for stacking. Dotted and solid lines show the InSAR pairs for ALOS and ALOS-2; red solid lines show those used for stacking in (c) in A1 and A2 are almost missing due to the averaged low-coherence area (Additional file 1: Fig. S1a). Moreover, we identified four robust subsidence areas (red arrows in Fig. 4a), corresponding to deforested agricultural fields (Fig. 2b), which usually experience thermokarst subsidence over a long period due to the melting of massive ground ice in this region. The magnitude of the subsidence signals for the four areas ranges from 0.5 to 3 cm yr −1 . Figure 4b shows the stacked vertical displacement rate using all selected SM3 interferograms (Table 1 and Fig. 4d). The signals in A1 and A2 are shown as subsidence with a rate of 1-2 cm yr −1 , but there is low coherence in A1 and A2 (Additional file 1: Fig. S1b). Moreover, we identified many subsidence areas around A1 and A2 at a rate of 0.5-1.5 cm yr −1 (Fig. 4b).

ALOS-2 SM3 stacking in 2015-2018
To derive the inter-annual surface displacement excluding seasonal components, we applied the stacking processing to four selected SM3 interferograms (Nos. 10-12 and 15 in Table 1 and Fig. 4d). Stacked interferograms using InSAR pairs that spanned the same seasonal interval but in different years were generated to minimize the effects of seasonal changes, such as thaw subsidence and soil moisture (Iwahana et al. 2016a). Thus, we could treat the signals in this stacking interferogram solely as inter-annual displacement. Figure 4c shows inter-annual subsidence at a rate of 0.5-2 cm yr −1 (T1-T5). Little displacement was detected in A1 and A2 in Fig. 4c.

ALOS-2 SM1 stacking in 2018
To see whether seasonal thaw subsidence was detected and whether more spatially small signals than that in SM3 were identified, we processed SM1 data obtained only in summer and autumn 2018 with short revisit intervals (14-42 days). Figure 5a shows the stacked Fig. 5 a Stacked vertical displacement between August 17 and September 28, 2018, derived from ALOS-2 SM1 data. The red circle shows smaller spatial-scale subsidence than that in Fig. 4. b Vertical displacement between August 17 and September 28, 2018 derived from a single InSAR (as in Fig. 3c). c Temporal distribution of InSAR pairs for stacking vertical displacement from August 17 to September 28, 2018 (Table 1 and Fig. 5c). The stacked interferogram displays much finer detail in displacement variations than those from FBD and SM3 because of the high spatial resolution of SM1 data. The magnitude of subsidence ranges from 0.5 to 4 cm during the period. Moreover, it should be noted that there are many subsidence signals on a spatially smaller scale (30 × 30 m 2 ) in the red ellipse. The stacked interferograms in Fig. 5a and the single interferogram of Fig. 5b (same as Fig. 3c) show similar spatial variation in surface displacement although the process varies for the same period (stacking vs. single pair). This indicates that this single SM1 interferogram (Fig. 5b) is high quality (i.e., produces little noise) and detects real subsidence signals similar to the stacking result, in which stacking is considered as a noise-reduced process (Fig. 5a).

Field observations in 2017-2018
To validate our results from the InSAR data, we performed leveling surveys in 2017 and 2018 at five plots (Fig. 2a, d). The results are summarized in Table 2. While there is little surface displacement at Plots A and D with a rate of − 0.2 ± 0.4 and 0.5 ± 0.6 cm yr −1 , Plots B and C showed significant subsidence with rates of − 2.6 ± 0.8 and − 3.2 ± 0.8 cm yr −1 . At Plot E, we obtained an interannual subsidence rate of − 0.5 ± 0.2 cm yr −1 during 2015-2018 as a minimum value.

Error sources on InSAR: interpretation and analysis
Although we consider seasonal and inter-annual displacements detected in this study to be caused by thaw settlement/frost heave including thermokarst, there are potential error sources related to InSAR phase bias caused by turbulent atmospheric effect and soil moisture change. However, the signals we identified in the stacking processes are considered as a more robust indication of true subsidence; this is described further in the following discussions.
The turbulent atmospheric effect may affect the interferograms. However, such signals could appear regardless of land cover and are unlikely to occur at specific open areas at the same time in the interferograms. With regards to localized atmospheric phenomenon in specific land covers in this region, it is known that mist occurs on summer days after sunset over open areas. The condensation of water vapor is caused by a difference in the heat capacities of forested and open ground. However, the height of the mist is usually less than 20 m, and the effect of the propagation delay of the SAR microwave due to the localized vapor is considered negligible. Moreover, similar mist occurs at almost all open (deforested) areas in this region and is unlikely to occur systematically at particular open areas. The effect by the mist may be due to the phase change of the constituent water from vapor to liquid droplets, which is not be due to the "turbulent" atmospheric effect. In either case, we confirmed each interferogram used in the stacking, and all signals in the discussed areas showed subsidence in multiple interferograms. Therefore, it is unlikely that these signals were due to the atmospheric effects.
With regard to the change in soil moisture and its related error, there is a possibility of apparent InSAR phase change in our results. The major causes of the possible errors have been considered to be the change in penetration depth of the microwaves, and the dielectric change that induces a wavenumber shift in the microwave propagation due to changes in soil moisture (e.g., De Zan et al. 2014; Zwieback et al. 2015). Zwieback et al. (2015) examined the dependency of the InSAR phase on soil moisture using airborne L-band SAR data, and concluded that the interferogram phase could change π/2 for a change in surface moisture of 20%. This change corresponds to 2-3-cm subsidence upon the 20% moisture increase at the L-band. In our study region, fluctuation in surface soil moisture is small because of the dry continental climate and the deep active layer. During the thawing season, as the thaw depth increases water holding capacity in the active layer will increase over time, and soil water in the near-surface layer can percolate into the deeper soil layers. The SAR data used in this study were acquired in the second half of the thawing season, and rainwater can easily infiltrate into the ground, and percolate down to deeper, dry soil layers in this region. In addition, the amount of summer precipitation in this area is very low (about 150 mm) due to the continental climate. In situ soil moisture data measured at inter-alas areas in late August and September from 2005 to 2018 indicate that the average soil water content of 10-cm thick surface layer was ~ 30%, and the amplitude of the fluctuation was ~ 7%. The InSAR phase error related to the change in soil moisture was estimated to be less than 1 cm. Fig. 2d, and E in Fig. 2a The value at A-D is a displacement rate for 1-year (2017-2018), and that at E for 3 years (2015)(2016)(2017)(2018)

Point
Measured values within 2σ (cm yr −1 ) In addition, error analysis of the stacking treatment in Rouyet et al. (2019) was conducted in our study. Rouyet et al. (2019) estimated the standard deviation of the stacking results of 0.25-0.35 cm per summer, assuming a standard deviation of 0.5 cm per interferogram due to the atmosphere, and using Eq. 11, given by Emardson et al. (2003). We carried out a similar analysis using Eq. 10 from Emardson et al. (2003), and the standard deviations of all interferograms. The standard deviation of the phase at each interferogram was calculated by masking out the displacements (indicated by the red arrows in Fig. 4a, and solid lines in Fig. 4b, c), which we identified using qualitative interpretation. Finally, the standard deviations of the stacking results of FBD (Fig. 4a), SM3 (Fig. 4b), and SM3inter-annual (Fig. 4c) were estimated to be 0.24, 0.11, and 0.23 cm yr −1 , respectively.

Spatial distribution of subsidence and land cover
The ALOS and ALOS-2 results reveal the spatial distribution of surface subsidence in Mayya (Figs. 3,4 and 5). The spatial variation of the displacement signals in FBD (Figs. 3a and 4a) is less clear than in SM3 (Figs. 3b and 4b,c). This is because the spatial resolution of SM3 (10 m) is twice that of FBD (20 m), and the B-perp of each SM3 interferogram was much smaller than that of FBD (Table 1; Ohki et al. 2018), leading to higher coherence (Additional file 1: Fig. S1a, b). The magnitude and location of the subsidence signals in the FBD and SM3 stacking results (Fig. 4a, b) differ for each period, which may indicate temporal changes in the subsidence rates. Moreover, the spatial variation of the displacement signals in SM1 (Figs. 3c and 6) is much clearer than that of SM3, owing to a shorter temporal baseline (14-42 days), and the SM1 images (3 m) having a much finer resolution than those of SM3. This enables us to identify the spatially smaller scale signals (~ 30 × 30 m 2 ) not detected by the SM3 images.
We compared the distribution of the displacement signals with land cover (Fig. 2b). The distribution of the subsidence signals corresponds to that for bare ground and/or grass, and not for forested areas. This agreement implies that the SAR microwaves can reach the surface of the ground and grass in open areas, but that it may not reach the surface in the forest due to interference of trunks and branches, which may mask ground displacement. When the tree structure deforms with or without the ground-surface displacement, those areas will show decorrelation with InSAR and information on ground displacement will not be available. Another possibility is that the forested ground was stable. The forest's vegetation layer can act as insulation to prevent permafrost from thawing (Iwahana et al. 2005;Shur et al. 2011), which may cause little surface displacement in the forested ground. The coherence in the forest was at least moderate (Additional file 1: Figs. S1 and S2), which indicates some information on the ground surface deformation is available in the interferograms. Therefore, we believe that some of the microwaves could reach the ground surface and return as coherent signals, providing information regarding ground surface deformation.

Site visit and implications
To confirm whether the InSAR displacement signals are valid, we visited some places where the subsidence signals were detected at the end of September 2018. Figure 6a is a local photograph at A1 in Fig. 3, which is located in an alas. Alasses are considered the final geomorphological stage of old thermokarst development (e.g., van Everdingen 2005), where no further thermokarst subsidence is expected. Thus, the signal at A1 may not be caused by the Fig. 6 Photographs at a the bottom area of the alas at A1 and, b Plot C in Fig. 2d. Polygonal texture is present at both the bottom area of the alas, and Plot C thermokarst. The two subsidence areas (A1 and A2) are shown in Fig. 3 (a single interferogram showing seasonal and 1-year displacement) and Fig. 4b, c (seasonal and inter-annual displacement derived from SM3) but are almost missing in Fig. 4a (seasonal and inter-annual displacement derived from FBD). During the ALOS operation, the alasses at A1 and A2 were almost flooded, due to meteorological conditions in 2006(Iijima et al. 2010. This led to a substantial decrease in coherence at A1 and A2 in most ALOS InSAR pairs, and the missing displacement value (Fig. 4a). In contrast, during the operation period of ALOS-2 (2015-2018), the ground surface in the alasses at A1 and A2 was mostly dry; therefore, ALOS-2 InSAR obtained high coherence and detected the displacements at A1 and A2 (Fig. 4b, c). In addition, the two subsidence signals in the stacking results appear clearly in Fig. 4a, b, but less so in Fig. 4c. This suggests that the signals in A1 and A2 are caused mainly by seasonal displacement. Figures 3a, c, and 5 show the seasonal subsidence at A1 and A2 of up to 4 cm from July to September in 2009 and from mid-August to the end of September in 2018, respectively. The signals in A1 and A2 in Fig. 3b indicate positive (i.e., uplift) because we considered the magnitude of uplift due to frost heave in 2017-2018 to be greater than that of the subsidence caused by a seasonal permafrost thaw in 2018. The result from Fig. 4c was derived using the interferograms including the result in Fig. 3b, but other SM3 interferograms indicated clear subsidence, which resulted in little displacement in A1 and A2 between 2015 and 2018.
It is interesting to note that we observed new developments of polygonal subsidence in the bottom area of the alas (Fig. 6a), which implies that massive ground ice remains in the permafrost under the alasses. Polygonal texture often emerges due to preferential ground subsidence along with the location of massive ground ice after surface disturbances and following thermokarst processes (Iwahana et al. 2016b). Although no inter-annual displacement is apparent to date (Fig. 4c), the land has the potential for further thermokarst progress (Ulrich et al. 2017).

Comparison of field observations and InSAR
We compared the results of field observations in 2017-2018 (Table 2, Fig. 6b) with the InSAR results (Fig. 7). Figure 7a, b shows an enlarged view of the SM3 single interferogram obtained on September 29, 2017 and September 28, 2018 (same as Fig. 3b) in the field survey area, and on October 2, 2015 and September 28, 2018 at Plot E. A comparison of the field survey and the InSAR-based displacement is shown in Fig. 7c. The three measured values (Plots A, D, and E) are comparable to those of the InSAR measurement (Fig. 7a, b), while the two measured values (Plots B and C) are noticeably greater than those measured by InSAR (Fig. 7a) by one order of magnitude.
Frost jacking of benchmarks could be a significant source of error in permafrost regions. In our study, detailed information about the installation of the reference pipe placed in the underlying permafrost for ground temperature measurement is not available. A potential Fig. 7 Comparison of surface displacement data obtained by field survey and InSAR. Enlarged views of SM3 single interferograms for a September 29, 2017, andSeptember 28, 2018 (as in Fig. 3b) in the ground survey area at Plots A-D, and b October 2, 2015 and September 28, 2018 at Plot E. c Comparison of the field survey and the InSAR displacement. The number of leveling measurements for Plots A-E are 33, 38, 36, 31, and 7, respectively error arises from possible frost jacking of the reference pipes, and resulting in the overestimation of the measured subsidence rate. However, the differences in surface displacements between Plots A and B, and Plots C and D are significant because the same reference pipes were used for both site pairs. Plots B and C were selected from areas of relative depression with bumpy relief, which is often observed as an initial stage of thermokarst development in this region (Soloviev 1973;Bosikov 1991). In contrast, Plots A and D were selected from relatively flat and stable areas surrounding Plots B and C, respectively. Judging by the average displacement of less than 1 cm at Plots A and D, the impact of frost jacking was negligible during our study, and the greater subsidence observed at Plots B and C was relative to the surrounding stable areas.
We also do not consider that the significant differences between the InSAR and field measurements to have been caused by turbulent atmospheric effects and soil moisture changes. Significant spatial variations in atmospheric moisture content within 100 m (the distance between Plots A and B or C and D) are unlikely. However, preferential changes in soil moisture at Plots B and C could have influenced the InSAR signals. We believe that the underestimate of subsidence by InSAR for Plots B and C could be partially attributed to a local decrease in soil moisture. The decrease in surface moisture at only Plots B and C is improbable considering the relatively depressed (concaved) relief.
Based on these discussions, there are three possible reasons for the discrepancy. Firstly, the InSAR measurements in this study might not have been sensitive to the spatial variations in surface displacement measured by our ground survey at Plots B and C Thermokarst is present locally within the observation area (Figs. 2e and 6b) and causes polygonal surface subsidence. Each polygon has a diameter of ~ 5-6 m (Fig. 6b), and the magnitude of thermokarst subsidence varies at the trough (the space between polygons) and center of the polygon, where the maximum and minimum subsidence occur, respectively. Given a spatial distribution of polygons in a pixel, most parts of the pixel are occupied by the center of polygons because the area of the trough is much smaller than that of the center. InSAR may not detect the subsidence in a trough but may detect averaged subsidence in the center of polygons, which may cause underestimation of surface subsidence by InSAR in thermokarst-affected areas. Second, field measurements might have been inaccurate, and third, the spatial representativeness of surface displacement by leveling does not match that of InSAR because of the insufficient number of measurement points at Plots B and C. However, it is unlikely that only the measurements at Plots B and C were inaccurate. The validity of the distribution and number of measurement points by leveling will be discussed in future work.
The inter-annual subsidence signals detected by SM3 (T1-T5 in Fig. 4c) cover approximately 400 × 400 m 2 , and we identified abundant polygonal relief in each area from high-resolution optical images. The T5 subsidence signal corresponds to Plot E, and the polygonal relief was also identified from the ground survey. The results of the ground survey are in good agreement with those of the InSAR (Fig. 7b, c) and the stacking (Fig. 4c). Assuming a similar situation, the subsidence signals should be significant for the T1-T4 areas although we have no field observations.

Separating seasonal and inter-annual changes in surface subsidence
Surface displacement related to the permafrost process is composed of seasonal (thaw settlement/frost heaving) and inter-annual (thermokarst) changes. Stacking processing is a simple method to extract small displacements assuming that the displacement rate is constant for each period. Therefore, our stacking results contain both seasonal and inter-annual displacements (Fig. 4a, b).
Major InSAR time-series analysis methods, such as the small baseline subset (SBAS) approach (Berardino et al. 2002), persistent scatterer InSAR (PS-InSAR) technique (Ferretti et al. 2001), and a combination of the two (Hooper 2008), were originally applied to derive small displacements such as inter-seismic deformation (e.g., Takada et al. 2018) and land subsidence (e.g., Ishitsuka et al. 2014). These methods are useful for identifying not only linear trends of surface displacement but also cyclic trends (seasonal changes) using fitting trigonometric functions and inversion algorithms. For example, Liu et al. (2015) used the SBAS technique to examine permafrost thaw subsidence in Alaska and derived thermokarst-induced and seasonal subsidence using 18 ALOS interferograms, with the inversion algorithm described in Liu et al. (2012). The inversion algorithm requires a sufficient number of interferograms to estimate an appropriate trend; however, the number of high-quality interferograms is often limited in some areas. Chen et al. (2018b) first applied the PS-InSAR technique to ALOS/ PALSAR data to derive permafrost thaw subsidence in the Qinghai-Tibet Plateau with a rate ranging from 0.3 to 3 cm yr −1 using 17 scenes of ALOS. SAR images acquired during the snow-covered season were used in the study, which may include errors associated with snow accumulation. Moreover, to extract PS points over the permafrost area, a low threshold should be set for determining the PS points, which may lead to low-quality results. In our study area, there are 17 scenes of ALOS/PALSAR data, and 12 of them were acquired in the snow-cover season. Thus, we only used five scenes obtained during the snow-free season to avoid the influence of snow accumulation. In ALOS-2/PALSAR-2, there were one or two acquisitions over Mayya until 2018 for a year, and fewer than ten acquisitions since the ALOS-2 launch. Therefore, the ALOS and ALOS-2 InSAR data (Table 1) available for time-series analysis are limited. Considering InSAR pairs for stacking, however, we derived the interannual displacement (Fig. 4c). This resulted from L-band data maintaining coherence over 3 years (Additional file 1: Fig. S2). We were unable to capture the complete of thaw and freeze cycle using interferograms (Fig. 5), and larger amounts of ALOS-2 data would enable us to simultaneously estimate seasonal and long-term surface displacement.

Advantages and limitations of L-band SAR for monitoring permafrost land
In the 2010s, many studies of deformation related to permafrost using SAR data have been reported (e.g., Liu et al. 2010Liu et al. , 2015Short et al. 2011Short et al. , 2014Iwahana et al. 2016a;Antonova et al. 2018;Chen et al. 2018a, b;Strozzi et al. 2018), which highlighted the advantages and disadvantages of different bands (X, C, and L). X-band data such as TerraSAR-X have a high spatial resolution and a relatively short temporal revisit interval of 11 days, enabling us to generate a highly coherent interferogram. However, the temporal decrease in coherence is more rapid than that for C-and L-band data; therefore, it was used to examine only seasonal changes. Antonova et al. (2018) used Ter-raSAR-X SAR data to examine seasonal thaw settlement in Northern Siberia, and they detected seasonal thaw subsidence of up to 2 cm in 2013, derived from eight cumulative interferograms with moderate coherence. They demonstrated that X-band InSAR is unsuitable for monitoring inter-annual subsidence due to a quick coherence drop for even 22 days, which has already been confirmed (Short et al. 2011;Wang et al. 2017). Sentinel-1A and -1B have a C-band SAR sensor with a revisit interval of 6 or 12 days, and a wide observation swath (up to 250 km for the nominal mode), which provides us much more frequent interferometric images than ALOS-2. The coherence of Sentinel-1 interferograms spanning 48 days is moderate, and capable of capturing spatial details and identifying non-uniform seasonal displacement (Chen et al. 2018a). Although coherent interferograms spanning 1 year, using images obtained at the end of summer could be generated using Sentinel-1 C-band data in some cases , our study indicates that L-band InSAR maintains coherence over a few years (Additional file 1: Fig. S2), which could derive inter-annual subsidence (Fig. 4b), even when using a small number of interferograms, because of the low latency of ALOS-2. These coherent interferograms spanning 1-3 years are necessary to measure inter-annual subsidence, such as thermokarst and better understand the thermokarst process Shiklomanov et al. 2013).
The moderate-to-high observation frequency of SAR satellites (shorter than 1 month) is crucial for deriving both seasonal and inter-annual displacement in permafrost dynamics. Although the 14-day revisit interval of ALOS-2 is much shorter than that of ALOS (46 days), there are some areas (e.g., Siberia) where the observation frequency of ALOS-2 is quite low. A new L-band SAR satellite, ALOS-4, will be launched in the 2021 Japanese fiscal year and will provide an observation swath of up to 200 km in high-resolution mode (3 m); therefore, the observation frequency will increase (Motohka et al. 2018). Moreover, ALOS-4 will fly in the same orbit as ALOS-2, which enables us to perform interferometry between ALOS-2 and ALOS-4 data (Motohka et al. 2018). Therefore, we expect more frequent acquisitions of L-band SAR data for a longer period. The currently operating SAtélite Argentino de Observación Con Microondas (SAOCOM) by Comisión Nacional de Actividades Espaciales (CONAE), and the upcoming National Aeronautics and Space Administration (NASA) and Indian Space Research Organization (ISRO) L-band SAR NASA ISRO synthetic aperture radar (NISAR) data will also contribute to a better understanding of permafrost dynamics. This is especially significant for the monitoring of inter-annual displacements, such as thermokarst, providing that ground observation data are obtained for a more accurate estimation of thermokarst subsidence by InSAR in Central Yakutia.

Conclusion
In this study, L-band InSAR data by ALOS and ALOS-2 were used to investigate surface displacement due to thermokarst in the continuous permafrost zone in Mayya, Central Yakutia. We focused mainly on spatial variation and temporal change in the surface displacement of deforested areas and alasses, and interpreted the detected signals, supported in part by ground survey. We also discussed the capability of L-band InSAR for monitoring thermokarst. From our findings, we derived the following conclusions: 1. Each interferogram generated from ALOS and ALOS-2 data shows subsidence and uplift in many areas. We identified numerous areas with subsidence signals in deforested areas that had been cultivated for farming in the 1970s, and alasses that developed thousands of years ago. In particular, the SM1 interferogram showed finer spatial variations in surface displacement and detected spatially smaller signals (30 × 30 m 2 ).