Toward the Detection of Permafrost Using Land-Surface Temperature Mapping

Abstract

Permafrost is degrading under current warming conditions, disrupting infrastructure, releasing carbon from soils, and altering seasonal water availability. Therefore, it is important to quantitatively map the change in the extent and depth of permafrost. We used satellite images of land-surface temperature to recognize and map the zero curtain, i.e., the isothermal period of ground temperature during seasonal freeze and thaw, as a precursor for delineating permafrost boundaries from remotely sensed thermal-infrared data. The phase transition of moisture in the ground allows the zero curtain to occur when near-surface soil moisture thaws or freezes, and also when ice-rich permafrost thaws or freezes. We propose that mapping the zero curtain is a precursor to mapping permafrost at shallow depths. We used ASTER and a MODIS-Aqua daily afternoon land-surface temperature (LST) timeseries to recognize the zero curtain at the 1-km scale as a “proof of concept.” Our regional mapping of the zero curtain over an area around the 7000 m high volcano Ojos del Salado in Chile suggests that the zero curtain can be mapped over arid regions of the world. It also indicates that surface heterogeneity, snow cover, and cloud cover can hinder the effectiveness of our approach. To be of practical use in many areas, it may be helpful to reduce the topographic and compositional heterogeneity in order to increase the LST accuracy. The necessary finer spatial resolution to reduce these problems is provided by ASTER (90 m).

1. Introduction

Permafrost, defined as ground that stays below freezing for more than two years, is an integral part of the cryosphere that is predicted to rapidly degrade under current warming climatic conditions (e.g., [1,2]). A quarter of the land surface of the Earth is subject to permafrost-related processes [3] and 23–26% of the land area in the northern hemisphere is permafrost [4,5,6]. The loss of permafrost affects the regional water balance and changes landscapes and ecosystems in cold regions. Erosion of soil due to melting permafrost not only poses a significant uncertainty in infrastructure stability but may also accelerate the release of carbon from the permafrost into the atmosphere [1,7].

Traditional mapping and monitoring techniques of permafrost are accurate but are presently labor-intensive and rely on interpolating small numbers of locally measured data across large areas in order to study regional-scale processes. Remote sensing addresses this difficulty by making spatially dense measurements over vast regions. Previous remote sensing studies have inferred permafrost extent from the instantaneous zero isotherm at the land surface or from changes in the dielectric marking the phase transition of liquid water to ice [8,9,10].

We propose to use land-surface temperature (LST) measured from space-borne thermal-infrared (TIR) sensors to map patterns in the annual cycle of LST in order to ultimately detect freeze/thaw surface dynamics. Previous studies demonstrated that in autumn or spring, the LST of moist/frozen soils remains at ~0 °C while the phase transition between water and ice occurs. This period is a well-known phenomenon called the “zero curtain” in permafrost studies (e.g., [11]) and its presence can be used to distinguish seasonally freezing ground from freezing dry ground, in which there is no zero curtain. The premise of our approach is that if the zero curtain lasts for several days, its presence can potentially be recognized and mapped with time series of daily LST images.

Repetition of the zero curtain for two or more annual LST cycles can be used to identify candidate permafrost pixels, and seasonal asymmetry in the duration of the zero curtain can likely be used to distinguish seasonally frozen soil from permafrost. Putkonen [12] suggested from modeling, substantially confirmed by Yi et al. [13], that the zero curtain should be present over near-surface permafrost during freeze-up, but not during the thaw, and vice versa in seasonally frozen soil. Over deep permafrost (~130 cm), de Pablo et al. [14] reported zero curtains in both seasons but shorter at the end of the freezing season. Although the zero-curtain duration is affected by seasonal changes in soil moisture and by the onset date of snow cover [13], this asymmetry should allow further winnowing of candidate pixels.

This paper does not reach so far as to identify permafrost, but simply establishes that precursor maps of the zero curtain are feasible to make with satellite-borne TIR imagers. We tested our approach in an arid region of the Atacama Andes near the Ojos del Salado volcano, close to the border between Chile and Argentina, where previous studies had documented the presence of frozen ground and near-surface permafrost [15,16,17,18].

1.1. Zero Curtain in the Surface and Soil Temperatures

The occurrence of a zero curtain in the annual ground-temperature cycle is attributed to the phase change of moisture in the soil due to the release of excess latent heat during the freezing and thawing seasons (Figure 1). Because keeping soil near 0 °C for several days generally requires moisture in the soil, the detection of the zero curtain implies the presence of soil moisture. Soil moisture may be ephemeral or seasonal, but it may further imply the presence of buried ice-rich permafrost. Putkonen [12] proposed that, in addition to the presence of moisture in the soil, the thermal gradient in the soil must become 0 °C m⁻¹ above the base of the active layer (depth of the thaw) for the zero curtain to form. It follows that the timing of the zero curtain indicates the presence of permafrost. During the thaw season, a high thermal gradient between the warming surface and the cold permafrost effectively conducts heat through the soil without forming a zero curtain; during freeze-up, the cold ground surface and the cold permafrost create a low thermal gradient, ensuring a condition at some depth in the soil where the temperature is at or near 0 °C for a prolonged period (Figure 1B). In other words, seasonally frozen ground not overlying permafrost may exhibit a zero curtain during the thaw, whereas soil underlain by permafrost may exhibit a zero curtain only during freeze-up.

Yi et al. [13] observed that the duration of the zero curtain revealed information about the permafrost, specifically the depth of the active layer. The heat-flow plateau at ~0 °C shown schematically in Figure 1B (from Putkonen [12]) would persist longer over thicker active layers but would break down more rapidly over thin layers. Measuring and mapping the duration of the zero curtain, therefore, would supply important information about the weather conditions, their change with altitude and geographic position, and changes over time relevant to climate change.

Significant characteristics of the zero curtain may be measured by three parameters: (1) The starting date of the zero curtain; (2) the duration of the zero curtain; and (3) seasonal asymmetry in the duration. While the zero curtain itself indicates the presence of soil moisture, its starting date reveals the beginning of the annual thaw and/or freeze-up and its change over the years may indicate the changing nature of the climate and seasonal fluctuations of soil moisture. The duration of the zero curtain may respond to the thermal regime of the soil and the presence of underlying permafrost and the depth of the active layer.

The analysis of daily measurements of LST is less sensitive to measurement precision and cloud cover than a single measurement of LST. The occurrence of the zero curtain in LST, regardless of the seasons, strongly indicates the presence of moisture in the soil and provides an opportunity to distinguish between dry and wet soil and map at a regional scale. Therefore, such a remote-sensing approach may be appropriate for detecting and mapping ground ice at moderate (10²–10³ m) resolutions.

1.2. Remote Sensing and Numerical Modeling Studies of Permafrost

Most of the traditional permafrost studies rely on the recognition of landscape changes in the field and extrapolating over large regions (e.g., [19]). Such field expeditions can be expensive and prohibitive in remote regions. If the changes in the landscape are large enough to be seen with remotely sensed data, the permafrost-related landscapes and their changes can be mapped at a regional scale. For example, Nitze et al. [20] used visible/near-infrared Landsat imagery (30 m resolution) from 1999 to 2014 to quantify the expansion of thermokarst lakes and the increased slumping of frozen coastal areas due to melting of permafrost. Such studies are invaluable to document landscape changes in permafrost areas but are not sufficient to explain the changes in the thermal regime of the soil.

Microwave sensors are sensitive to soil moisture and temperature (e.g., [21]), and exploiting the fact that freezing changes the dielectric constant of moist soils, numerous studies of radar measurements have been used to map the state of freeze and thaw in the soils (e.g., [22,23]). Such regional-scale mapping approaches based on microwave sensors are useful to infer the changes in seasonally frozen grounds. However, the existence of permafrost underlying the seasonally frozen ground has not been inferred from just the dielectric constant or the state of freeze/thaw, and the resolution of the microwave sensors (~6–60 km) can be a limiting factor to detect small changes in lateral extent of permafrost areas. With the goal of mapping the depth of the active layer, Liu et al. [24] used interferometric synthetic aperture radar (InSAR) to monitor surface deformation over permafrost. Because ice and water have different densities, the ground surface settles during the thaw and the amount of surface subsidence can be used to estimate the thickness of the active layer if the vertical distribution of pore water within the soil is known. Wang et al. [25] used multi-temporal TerraSAR-X backscatter intensity and interferometric coherence along with the relationship between vegetation cover and the permafrost as the basis for classifying and mapping permafrost landscape features. Chen et al. [26] proposed a method to estimate the active layer thickness using time-series P-band polarimetric synthetic aperture radar (SAR) observations, achieving retrieval with errors <0.1 m in some cases. Yi et al. [13] combined field observations, active layer thickness, and soil moisture maps derived from low-frequency (L + P-band) airborne radar measurements, and global satellite environmental observations to investigate the sensitivity of the active layer to recent climate trends.

Thermal infrared (TIR) sensors provide direct measurements of skin temperatures of land surfaces. Widely used LST products are derived from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) onboard the Terra spacecraft, the Thematic Mapper (onboard Landsat), and from the Moderate Resolution Imaging Spectroradiometer (MODIS) onboard two satellites: Terra (with a morning overpass) and Aqua (with an afternoon overpass). Terra, Landsat, and Aqua share a common sun-synchronous orbit with different overpass times. Among the TIR imagers on these satellites, only MODIS provides daily measurements of LST. Numerical modeling studies based on the empirical relationship between the ground temperature and air temperature have taken advantage of the MODIS LST products (e.g., [27,28,29]) to derive static maps of the spatial distribution of permafrost. Recently, a numerical modeling study has combined the MODIS LST products with gridded air temperature data to derive time-series maps of permafrost extent [5]. These dynamic maps of changes in permafrost are important to assess the impact of a warming climate to the state of the permafrost. The approach proposed in this article is a contribution toward creating temporally and spatially dynamic maps of the quantified state of permafrost and seasonally frozen ground, complementary to the previous regional-scale studies.

2. Materials and Methods

2.1. Remotely Sensed and Climate Data

MODIS provides daily LST products (MOD11A1 and MYD11A1) at 1-km spatial resolution using TIR data measured from the Aqua and Terra platforms, respectively [30,31]. We used MYD11A1 (Aqua) daytime LST products from 2003 to 2017. The MODIS LST is estimated using the split-window algorithm [32] with an estimated uncertainty of ~4.5 °C [33]. We compiled the daily MODIS LST over the study area and converted the data for processing with MATLAB. Daily coverage is possible using oblique views from adjacent overpasses. The local time for MODIS measurements in the LST compilations, therefore, ranged from 1:10 to 2:50 PM.

We used an ASTER Surface Kinetic Temperature (AST_08; 90 m pixel⁻¹) and visible/near-infrared (VNIR) ASTER products (AST_L1B; 15 m pixel⁻¹) to assess the variability of the land surface and the related temperatures within a single MODIS pixel, and to check the distribution of snow. The ASTER products were derived from a scene acquired on 8 April, 21 August, 22 September, and 8 October of 2017 at 11:43 AM local time [34].

Air temperature variations over different altitudes were evaluated using the long-term (1981–2010) monthly mean air temperature data from the Global Historical Climatology Network and Climate Anomaly Monitoring System (GHCN CAMS; 0.5° pixel⁻¹ [35]). These low-resolution air temperature data were resampled at a 30 m pixel⁻¹ resolution using the SRTM 1 arc-second DEM [36] and the long-term (1981–2010) monthly mean lapse rate derived from the National Centers for Environmental Prediction and the National Center for Atmospheric Research reanalysis data (NCEP/NCAR; 2.5° pixel⁻¹ [37]). We followed the procedures described in [38] to resample the air temperature data to a higher spatial resolution to see the fine spatial variations of air temperature.

2.2. “Threshold Window” Filtering Algorithm

We developed an algorithm in the MATLAB environment that analyzes the daily MODIS LST during the first and second six months of the year at each pixel, and tests whether the daily LST satisfies the following rules indicating the occurrence of a zero curtain:

(1)

LST must be between −3.5 and 3.5 °C (defined below as the LST during a zero-curtain event, or “zero-curtain LST”);

(2)

The number of consecutive zero-curtain LSTs must be >3;

(3)

The number of days with missing data between two identified zero-curtain LSTs must be <3; and

(4)

The total number of zero-curtain LSTs must be >5.

In other words, we split the annual daily LST data into two halves, each containing the thaw or freeze-up period. In each half, we counted the number days during which the −3.5 < LST < 3.5 °C. We identified a zero curtain if the number of consecutive days with −3.5 < LST < 3.5 °C exceeded five. At the end of the counting iteration for a pixel, the starting and end date for the zero curtain were recorded, from which we calculated the difference, which we defined as the duration of the zero curtain. Currently, no correction is made for snow cover. Our algorithm returns as output for every pixel matrix for the zero-curtain duration and the starting dates for thawing and freezing seasons for each year. The output matrix was converted back to georeferenced raster files using the same gridding scheme and spatial resolution of MYD11A1 products. All maps and raster data produced from this study are available in georeferenced TIFF files, as well as the detailed descriptions of the data, in the online supplement to this article (Table S1).

2.3. Validation Sites

We chose three validation sites near the Ojos del Salado volcano in the Andes east of the Atacama Desert in Chile (Figure 2). At least one 1-km² area surrounding each study site is homogeneous in terms of composition and small-scale roughness, topographically flat, and sparsely vegetated. Climate data are given in

Table 1. Validation sites, their locations, and local climate. Precipitation data are too low resolution to show differences among the sites.

Site	Location	Altitude, m asl	Mean Annual Climate Parameters
			Precipitation 1, mm w.e.	Air Temperature 2, °C	Lapse Rate 3, °C km−1
#1	27°02′54” S, 69°04′52″ W	3815	58	−0.8	5.5
#2	26°57′31” S, 68°49′09″ W	4415	58	−3.1	5.5
#3	27°00′09” S, 68°42′55″ W	4910	58	−5.0	6.3

¹ All sites show the same precipitation due to 0.5° resolution of GPCC V7 data [42]; ² [35]; ³ [37].

. Only ~1% of the total surface area over the whole region including the validation sites is covered with shrubs, grassland, water bodies, and permanent snowfields or glaciers (supplementary Table S2). Two of the sites were chosen in areas that were previously found to have periglacial features [39] and possibly permafrost [16,40]. However, intermittent or thin snow cover on Ojos del Salado implies complicating factors due to increased soil moisture during the melt season and the possibility of near-surface ground ice. The sites were chosen at different altitudes, ranging from 3815 to 4910 m asl, to assess the variability in the zero curtain at various altitudes. The lowest site lacked periglacial features and was intended to act as control for the higher ones.

Site #1 is located on one of the low-gradient (1°) alluvial fans (~140 km²) emanating from the rivers Quebrada Cienaga and Rio Lomas (Figure 2). These rivers form braided channels that lead to Salar de Maricunga, but the channels remain dry most of the time. The surface of the fan is paved with rhyodacitic/andesitic (cf. Baker et al. [41]) gravel 2–3 cm in diameter and 1 cm thick, overlying silty sand. A few of the clasts are 3–5 cm in diameter and <2 cm thick. Site #2 is located on an alluvial fan (~13 km²) below the volcanoes Nevado Tres Cruces. The surface of this fan is a low-gradient (~2°) rhyodacitic/andesitic/dacitic gravel pavement similar to Site #1, with gravel 2 cm in diameter and 1 cm thick. Site #3 (~1 km²) is on a steeper alluvial fan (~10%) in the Valle de Barrancas Blancas and has the most complex surface in terms of composition (pyroclastic, rhyolitic, dacitic) and topography. It is covered with gravel 3–5 cm in diameter and <3 cm thick. Some large wind-polished boulders are scattered on the surface. Almost-parallel, shallow, dry channels ~20 cm deep dissect the surface. Small grasses and flowers may cover the site in summer. At all three sites, the percent of the surface covered by gravel was >90%. Photos showing the surface features at the validation sites are provided in supplementary Figures S1–S7.

Surface features in the study region indicate periglacial processes and seasonal freeze and thaw. For example, in the Valle de Barrancas Blancas (Figure 2), sorted and unsorted patterned ground, surface extrusion due to cryoturbation, rock glaciers, cryoplanation surfaces, and slopes affected by gelifluction have been observed by Buchroithner and Trombotto Liaudat [43]. Nagy et al. [16] measured ground temperatures at six sites at different altitudes, ranging from 4200 to 6890 m asl, around Ojos del Salado during 2012–2016 to show that the seasonal freezing front reaches depths of at least 35–60 cm at these altitudes. Zero curtains were observed in both spring and autumn for two more weeks at all six sites, except for the two highest sites, at 6750 and 6890 m asl, where the ground temperatures remained below 0 °C for most of the year (see Section 2.4 below). On the basis of these measured ground temperatures, Nagy et al. [16] concluded that permafrost likely exists above 5250 m asl. The likelihood of permafrost at these altitudes is consistent with suggestions based on the mapping of permafrost and periglacial landforms [39,44] and numerical modeling of air temperatures [40]. Nagy et al. [16] concluded that it is unlikely for permafrost to exist below 4550 m asl.

2.4. In-Situ Measurements at Validation Sites

In November 2016, we installed TMC6-HC temperature sensors with an accuracy of ±0.25 °C at each site to measure subsurface ground temperatures every hour at depths of 2, 10, 20, and 40 cm (Figure 2; Table 1). The data from the thermometers were collected using HOBO^® H-8 4-channel data loggers. The near-surface sensor was for comparison to the MODIS LST data, averaged over 1 km². The deeper subsurface sensors were to evaluate the 2-cm temperature in terms of the diurnal cycle. When we removed the probes (March 2018), we examined the soils. See soil descriptions at Site #1 (Table S3), at Site #2 (Table S4), and at Site #3 (Table S5) in the online Supplementary Materials.

3. Results

3.1. The Zero Curtain in the LST Data

The LST for the ~2 PM MODIS overflight of Ojos del Salado during the thawing season of 2017 is shown for Site #2 in Figure 3. The plateau near 0 °C between June 20 and July 14 is the zero curtain. 2 PM is near the peak of the diurnal temperature cycle, and thin ice films forming at night will likely not register then. The lack of LST records reflects the cloudy days.

MODIS LST images were consistently recorded for the three validation sites between 2003 and 2017. The seasonal temperature cycle is clearly visible in Figure 4, although the short zero curtains are not evident at this compressed temporal scale. The persistence of low winter temperatures was less at all three sites than the persistence of high summer temperatures. Gaps in the record are times when clouds obscured the ground. The maximum temperatures ranged among sites from ~50 °C at sites #1 and #2 to ~35 °C at Site #3. The minimum temperatures ranged much less, only rarely dipping below 0°C during the afternoon even at Site #3 (4910 m asl), 460 m above the altitude at permafrost was thought possible by Nagy et al. [16]. The LST at the lowest site, Site #1, consistently stayed above freezing for almost all years, consistent with the low probability of permafrost existence based on the surface mapping and numerical modeling of air temperatures [39,40,44]. The highest site, Site #3, near Barrancas Blancas shows consistent zero curtains during both freezing and thawing seasons.

3.2. Occurrence of Zero Curtains in the Study Area

We analyzed 3393 pixels (87 × 39 km grid) with a 1-km spatial resolution in the study area for each of the 15 years (2003–2017). Approximately 60% of the study area is located above 4550 m asl, the altitude above which permafrost likely exists per Nagy et al. [16]. We produced two types of maps for thawing and freezing seasons of the years 2003–2017 using the Threshold Window algorithm: (a) Number of days with zero curtain (duration), and (b) starting date of the zero curtain. The number of pixels with detected zero curtains (Figure 5) shows significant differences between the thawing and freezing seasons. The area with zero curtains during the freezing season remained largely constant throughout the years. However, the area with zero curtain during the thaw has increased over the years. A notable exception was during 2015, where the number of pixels with zero curtain during the freezing season was anomalously high.

Zero curtains were observed in ~5% to 40% of areas above 4550 m asl between 2003 and 2017. In the majority of cases, the zero curtain was observed only during a single season. In each year, except 2015, less than 50 pixels showed zero curtains in both seasons. In other words, there was very little overlap between the areas with zero curtain detected only during freeze-up and the zero curtain detected only in the thaw (supplementary Figure S8). During the thaw, the zero curtain was detected largely in areas at low altitudes (e.g., Figure 6).

During freeze-up between 2003 and 2017, zero curtains were detected in late April to late May (supplementary Figure S9). In the higher altitude areas, the zero curtain started earlier than at lower altitudes, suggesting that the MODIS LST measurements reflect the temperature decrease with increasing altitude. During the thaw, zero curtains were detected between late June and early October. Zero curtains during the thaw of 2003–2010 were detected mostly in August. However, after 2010, zero curtains were detected earlier, during late June to the middle of July (e.g., Figure 7). These earlier starting dates during the thaw seasons may reflect the recent warming trends in the air temperatures over the Andes [45].

3.3. Comparison Between MODIS LST and In Situ Measured Ground Temperatures

The MODIS daytime LST is a measure of the ~2 PM skin temperature integrated over each 1 km² pixel. The skin temperatures are subject to short-term fluctuations caused by gusts of wind and the shadows of high-altitude clouds located over adjacent pixels. The km-scale skin temperature is weighted by the temperatures of the gravel of the pavements on the fans, whereas the subsurface temperature probes record the temperature of the soil beneath the pavement. Despite such different representations of ground temperatures, the measurements of daily MODIS LST and subsurface ground temperatures at the 2 cm depth in the three validation sites during 2017 agree well during the warm seasons (Figure 8), demonstrating the robustness of MODIS LST measurements for monitoring purposes. However, the correlation is less well established during the cold seasons. During 2017, the ground temperatures at the three validations sites were consistent at all depths (2–40 cm) during May–September, suggesting that the surface was covered with snow (see MODIS snow cover data, MYD10A1 [46], in Figure 8). Both MODIS snow cover data and the in situ measurements of ground temperatures at Site #3 show the fewest number of days with snow cover.

The MODIS LST at sites #1 and 2 show 15 days of sustained zero curtain during the thawing season of 2017. However, during this time at both sites, the subsurface temperatures measured at 2 cm were ~5–10 °C lower than the LST values. Only at Site #2 did the ground temperatures at the 40 cm depth show a persistent zero curtain, during which time the MODIS LST appear to agree (Figure 8B). We discuss in Section 4 below the complexities in the TIR remote sensing and their implications for the MODIS LST data.

4. Discussion

4.1. The MODIS LST Product

The basis for the detection of the zero curtain in this study is the MODIS LST product, of low (1 km) spatial but high (daily) temporal resolution. The MODIS LST is constructed from MODIS bands 31 and 32 (10.78–11.28 and 11.77–12.27 µm) using a split-window algorithm that corrects the data for atmospheric effects. Calculating LST from atmospherically corrected surface-emitted radiance is sensitive to the thermal emissivity: An error of 0.01 at 11 μm (MODIS band 31) corresponds to a temperature error of 0.6 K in the inversion of the Planck equation for the single wavelength. The default surface type (and emissivity in MODIS band 31) is soil (0.97: [47]). Although the MODIS LST algorithm assumes emissivity values for different surface types [48,49,50], only two are relevant in our validation sites on gravel pavements: Rock and snow (band 31: ε = 0.98). Coincidentally, the dacitic and volcanic pebbles composing the surfaces of the validation sites have the same emissivity at 11 μm used by the MODIS LST algorithm for soil. Therefore, the composition of the geological surface is unlikely to introduce temperature errors in the output images. MODIS snow-cover products (MOD10 and MYD10) are integrated into MODIS LST so that the calculation uses the proper emissivities, whether for the gravel surface or snow. Any discrepancy in LST due solely to the emissivities of the geologic surface and snow cover would be too small to detect readily in the time histories of LST for the three validation sites (Figure 4).

The cloud cover is detected and masked out in making the MODIS LST products [51] so that only the unobscured ground surface is represented. However, the lack of daily LST data during the seasonal thaw and freezing seasons may hinder the opportunities to detect the zero curtain occurrence. The percentage of data loss due to clouds over the study region from 2003 and 2017 ranged from ~20% to 60% in the mountainous areas. The percentage of data loss during March–September from 2003 to 2017 at the validation Site #1 ranged from 16% to 42% whereas at sites #2 and 3 it ranged from 27% to 56% (supplementary Table S6 and Figure S10).

The MODIS LST represents the dynamic balance between heating due to absorption of sunlight, diffusion of heat into the soil, and cooling due to emitted thermal-infrared radiance. The ‘true’ kinetic surface temperature, on the other hand, is susceptible to additional environmental changes, such as wind, evaporation, or condensation of water. These typically will make the LST more variable than the subsurface measurements, which are buffered by a large mass of adjacent soil. During the zero curtain times in 2017, the difference between the MODIS LST and the in situ ground temperatures at the 2 cm depth ranged between 1 and 6 °C at the validation sites (Figure 8). Considering the accuracy of the MODIS LST product of ~1 K [52] and the in situ ground temperatures, especially at the 2 cm depth, were comparable for much of the year. Exceptions to this rule evident in Figure 8 are discussed below.

4.2. MODIS LST and Subsurface Temperature Profiles

Figure 8 compares the LST to the subsurface kinetic temperatures measured at the same time at the three validation sites. Although the LST and subsurface soil temperatures differ significantly depending on the type of vegetation covers (e.g., [22]), the surface at the validation sites remain barren most of the time (e.g., Table S2). As expected, the temperature at the 2 cm depth generally tracks the LST closely at all three sites, although it is lower by up to ~10 °C in summer. This may be due to the phase lag as the solar heating wave slowly penetrates the ground. In winter, the 2-cm temperature lowers to match the other subsurface values, but the afternoon LST rarely drops below zero.

Snow is highly reflective in the visible part of the spectrum (50–80%), yet 5% of the incoming solar energy is transmitted through cover 10–20 cm thick [53]. In contrast, snow absorbs efficiently in the thermal-infrared spectrum, such that the emitted radiation may originate only in the top 2 mm of the snow surfaces [54,55]. The nearly isothermal character of the subsurface as shown in Figure 8A is likely due to an insulating blanket of snow 10 cm or more in thickness, allowing the shallow subsurface to lose heat to deeper soil even as the snow surface is heated by the sun. In the limit, the LST of the snow cover could exceed the freezing point by at most a few °C, yet temperatures as high as 34 °C were recorded, even as the subsurface temperatures continued to decrease. This suggests that numerous large, sun-heated cobbles protruded from the thinning snowpack or that unresolved patches of snow-free ground had appeared but not over the subsurface sensors. By mid-October, the 2-cm temperature values have risen to match the LST, likely as the snow melted.

At Site #2 (Figure 8B) from mid-July until late September, the LST was lower than the temperature 2 cm down. The MODIS snow product indicates that there was snow cover during this time. This situation is the opposite of what occurred at the same season in Site #1. How could this happen, and what does it imply for the detection by MODIS of a zero curtain?

A likely explanation is that the snow cover was a thin veneer, thick enough so that the thermal infrared signal was radiated from its surface, unaffected by the ground beneath. On the other hand, the veneer was thin enough (i.e., <10–20 cm) that the ground was heated by irradiance from the sun at visible wavelengths. Thus, in mid-May, both the LST and the 2 cm temperature lowered to a few degrees below zero as the first snowfall was recorded.

During the time the ground was snow-covered or obscured by cloud, no meaningful zero curtain could be detected. This period included the likely freeze-up. By early June, the cloudiness cleared, and the surface of the snow rose from −10 °C to zero as the sun heated and thinned the snow cover. From mid-June to mid-July, the LST indicated a zero curtain. Thus, misleading zero curtains that pertain to thawing snow can be returned by the “Threshold Window” filtering algorithm.

At the high-altitude Site #3 (Figure 8C), subsurface zero curtains (10, 20, and 40 cm depth) occurred during both thaw (early May) and freeze-up (early October), but none occurred at the 2 cm depth or at the surface, both of which were 20 °C or warmer. Furthermore, except for the temperatures at 2 cm, the subsurface sites had similar temperatures throughout the winter, between the zero curtains. MODIS reported no snow cover before June or after August, so that snow has no role in explaining these relationships then. During the winter, however, on individual days, temperatures at both 2 cm and the immediate surface dropped ~12 °C to near zero, an occurrence likely associated with ephemeral snow cover. During the summer, temperatures at 20 and 40 cm remained close to each other while temperatures at 10 cm rose to intermediate values between the temperatures at the deeper sensors and the 2 cm depth. These observations likely resulted from solar warming of the surface. At Site #3, there was sufficient soil moisture to allow zero curtains at depth, but the circumstances precluded any at the surface.

4.3. Effects of Scene Roughness

The alluvial fan surfaces at the validation sites were compositionally homogeneous and flat on a coarse km scale. However, we noticed that local fine-scale topography was variable, ranging from dry channels to small periglacial patterned grounds (supplementary Figures S1–S5). The manifestation of topographical variation in the emitted radiance at the thermal infrared spectrum can be seen in Figure 9, where subtle topographic features not discernible in the photograph can be resolved as bright pixels in the FLIR image. Inclusion of these features in low-resolution pixels causes the histogram of temperatures within those pixels to broaden. For example, the standard deviations within pixel aggregations simulating lower-resolution data in the plain in Figure 9 range from about 0.2 K at 2 m pixel⁻¹ to 2.5 K at 380 m pixel⁻¹. Topographic effects influencing the effective surface temperatures at the ASTER or MODIS scales are not large enough to affect the analysis in this study significantly. However, during the thaw, the radiance from unresolved patches of snow will mix with that from adjacent bare soil to lower the effective daytime LST.

4.4. Spatial Resolution and Radiance Mixing

Although the validation sites appear to be homogeneous from the standpoint of the MODIS LST data, analysis of zero curtains elsewhere will encounter more complex surfaces and in any case there may be seasonal patches of snow unresolved by MODIS. ASTER, at 90 m pixel⁻¹, will produce more unmixed (“pure”) pixels in complex terrain, and thus will produce more pixels in which reliable assessments of the zero curtain can be made (Figure 10). Small seasonal snow patches provide one example, important in this study, of temperature mixing that may be reduced with improved resolution.

For this study, the temperature window in the Threshold Window filtering algorithm was set at ±3.5 °C to account for the MODIS LST accuracy of ~1 °C (1 σ), any variations in emissivity of different types of rocks, and effects of local topographic features that cannot be resolved at the spatial resolution of MODIS LST data. Many but not all mixed pixels of soil and snow will be filtered out with this setting.

Adjusting the size of the Threshold Window decreases or increases the number of data points to be counted as part of a zero curtain, decreasing or increasing false negatives or false positives, but a better approach might be to reduce the pixel size. For example, horizontal changes in the extent of permafrost in Tibet Plateau can be estimated to ~460–920 m in 30 years after accounting for the maximum vertical change of permafrost base during the same period (80 m [56,57]) and the average slope of the interior of the plateau (~5–10° [58]). The threshold window algorithm cannot capture changes at such a 100 m scale due to the relatively low resolution of the MODIS LST data. This limitation might be overcome by using ASTER data, with its higher spatial resolution than MODIS, but repeat acquisitions of ASTER images are too infrequent. A compromise may be possible with data fusion: Surface compositions can be mapped using the 15-m pixel−1 visible and 30-m pixel−1 near-infrared data while variations in the surface emissivity can be resolved from the 90-m pixel−1 ASTER products (Figure 11). With this information and the daily LST, it may be possible to generate the higher 90-m spatial resolution of the ASTER surface temperature product using the STARFM algorithm [59]. Data fusion might also be used with the high spatial but irregular and temporal resolution data of ECOSTRESS (

Table 2. Comparison of thermal-infrared sensors.

Sensor	Resolution		Launch Date
	Spatial	Temporal
GOES	Low (4 km)	High (3 h)	1981
AVHRR (NOAA)	Moderate (1 km)	High (daily)	1979
MODIS	Moderate (1 km)	High (daily)	2000
ASTER	High (90 m)	Low (16 days)	2000
Thematic Mapper	High (60 m)	Low (16 days)	1982
ECOSTRESS 1	High (38 × 69 m)	1 h of science data day−1	2019

¹ The ECOSTRESS products are in development and data are available only through “Early Adopters” program.

). However, to monitor subtle changes in permafrost extent in the future, we likely will need data with both high spatial and high temporal resolution.

4.5. Seasonal Zero-Curtain Duration

Putkonen [12] stated that there is an asymmetry in the occurrence of the zero curtain during spring and autumn over permafrost, where the zero curtain appears during freeze-up, and over seasonally frozen ground, where the zero curtain appears during thaw. Figure 12 examines this asymmetry both as a correlation between the duration of freeze-up and thaw and as a function of altitude in the validation sites during 2015. Above ~5200 m asl, the zero-curtain duration of both freeze-up and thaw is about the same; at lower altitude, the duration of freeze-up exceeds the duration of thaw~40%.

Is the weak asymmetry revealed in Figure 12A due to a systematic trend, or is it the result of weather anomalies? In late March of 2015, before the general freeze-up, an extreme heat event occurred in central and northern Chile followed by an anomalously high precipitation event (>20 mm day⁻¹ [60]). Compared to other years, this extreme weather event must have supplied enough moisture to the soil to have resulted in the anomalously high number of pixels with zero curtains during both freeze-up and thaw that year (Figure 12A), so the anomaly did not increase any apparent asymmetry. During freeze-up, the zero curtain was detected in early April in high-altitude areas and in early June in low-altitude areas. The zero curtain during this time persisted for about two weeks. Later in the year, during the thaw, the zero curtain was detected only around early July with no contrast between the high- and low-altitude areas. The anomalously high precipitation from the March event may have been preserved in the soil over the winter and the increased soil moisture may have contributed to the sustained zero curtain during the thaw of 2015.

The weather anomaly aside, Figure 12B shows that the in the study area the maximum duration of the zero curtain during freeze-up, but not during thaw, may have increased erratically over the study period. The increase during freeze-up suggests that changes in the zero curtain can be detected on the decadal scale. The more constant and smaller number of zero curtain days during thaw indicates that the results are not random and that there has been little change. Taken together, the data plotted in Figure 12 suggest that there is little evidence at Ojos del Salado for the strong asymmetry predicted by Putkonen [12] for the zero-curtain freeze-up/thaw duration over permafrost. It may be that summer, but not winter, precipitation has been increasing in recent years but is not sufficient to carry over to the thaw. In all, these results are encouraging for the use of the zero curtain in climate-change studies.

4.6. Potential for Quantified Mapping of Seasonally Frozen Ground and Permafrost

The occurrence of the zero curtain requires moisture in the soil and low temperatures. Therefore, the identification of a zero curtain is a useful indicator of ice-rich frozen ground and possibly permafrost. We identify the following four parameters for the zero curtain that can be quantified using the Threshold Window algorithm: (a) Starting date, (b) duration, (c) seasonality, and (d) persistence over more than two seasons. The changes in the starting date of a zero curtain over time may indicate changes in the thermal regime of the soil due to changing environmental and climate conditions. According to Putkonen [12], the thermal regime in the soil would be different depending on the existence of permafrost underneath, and the asymmetric seasonality of the zero curtain may indicate whether it is due to seasonally frozen ground or the possible existence of permafrost. The intensity of freezing or thawing may dictate the duration of zero curtain, the changes of which can be mapped over time. On the basis of these parameters and the criteria for classifications of permafrost conditions (e.g., [61]), it is possible that a quantitative map of permafrost extent and the depth of the active layer can be estimated.

The probabilistic horizontal extent of permafrost based on the empirical relationship between ground temperature and mean annual air temperature has been mapped (e.g., [40]). The detection of permafrost using the Threshold Window algorithm, therefore, can validated using the probabilistic maps of permafrost produced from climate and geomorphological data. Further, the depth of the active layer can be estimated by numerically modeling heat diffusion into the ground using surface temperatures estimated from ASTER or MODIS data as constraints. For example, a heat-transfer model of soils that accounts for different estimates of moisture content [62,63] was run using MODIS LST (8-day mean) as an initial condition and used to numerically construct the thermal profiles of soils in a well-known permafrost region of Tibet (Figure 13; [64]). The ground temperatures calculated at various depths over time allowed the depth to the zero isotherm to be estimated for 2002 and 2009, consistent with the depths to the base of the active layer measured in the field [65]. Figure 13 shows that the depth of the zero isotherm decreased about 20 cm over the 7-year period. The model suggests that there was an extended period when −5 °C < LST <5 °C that may have been a zero curtain (intersection of the dark blue color and the top of the plot) in April–May 2002, and that it was shorter in 2009. Similarly, a zero curtain was predicted for October, and it too was shorter in 2009. Thus, this model suggests that there may have been a zero curtain not just at freeze-up but also during the spring thaw, inconsistent with Putkonen [12] (Figure 4). Further testing of the model will be necessary to resolve this issue.

5. Summary and Conclusions

We proposed a new approach to detect the occurrence of the zero curtain in land-surface temperature images, maintenance of the temperature at or near 0 °C due to the release of latent heat from freezing or thawing moisture in the soil. We used a “Threshold Window” filtering algorithm to analyze daily MODIS land surface temperature (LST) data over a cold region of the Atacama Andes in Chile, where periglacial features have been documented. Our results demonstrated that the zero curtain could be consistently identified in the MODIS data. This demonstration opens the path to identifying ice-rich permafrost using time series of LST images.

The duration, seasonality, and the starting time of the zero curtain for each year between 2003 and 2017 were mapped. We measured subsurface ground temperatures at depths of 2–40 cm at three validation sites located at altitudes of 3815, 4415, and 4910 m asl in 2017. The in situ observations of subsurface temperatures complemented the MODIS LST records and showed that, with exceptions during the spring, the LST agreed within a few degrees with temperature at the 2 cm depth. The subsurface temperature records showed clearly the presence of the zero curtain at depth, out of phase with the zero curtain at the surface, and illustrated the effects of thinning snow cover on the recovered zero curtain. Over the period of the study, the duration of the zero curtain during the thaw and freeze-up seasons appeared to fluctuate, consistent with what might be expected from annual changes in weather conditions. However, on Ojos del Salado, the zero curtains at freeze-up tended to be longer than those at thaw, and also may have shown a tendency to lengthen over the study period. We did not find convincing evidence of permafrost in our data or in the field, although there was local evidence of segregation ice.

Our test sites were chosen to minimize vegetation cover that would interfere with satellite measurement of the land-surface temperature, but in practical applications, vegetation cover and snow may hinder the availability of remotely sensed LST for this approach. Small-scale variations in surface composition and topographic conditions create different thermal components that are integrated over large areas for products with moderate spatial resolution, such as MODIS LST. Thermal-infrared sensors with higher spatial resolution, such as ASTER and ECOSTRESS, can be used to resolve the 100-m-scale horizontal changes in the extent of seasonally frozen ground and permafrost. Numerical modeling of the thermal regime of the ground using these surface temperature records may be useful to estimate the depth of thaw in ice-rich frozen grounds and permafrost.