Next Article in Journal
A Modular Processing Chain for Automated Flood Monitoring from Multi-Spectral Satellite Data
Next Article in Special Issue
Evolution of Near-Shore Outwash Fans and Permafrost Spreading Under Their Surface: A Case Study from Svalbard
Previous Article in Journal
Integration of ZiYuan-3 Multispectral and Stereo Data for Modeling Aboveground Biomass of Larch Plantations in North China
Previous Article in Special Issue
Seasonal Progression of Ground Displacement Identified with Satellite Radar Interferometry and the Impact of Unusually Warm Conditions on Permafrost at the Yamal Peninsula in 2016
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Assessing Spatiotemporal Variations of Landsat Land Surface Temperature and Multispectral Indices in the Arctic Mackenzie Delta Region between 1985 and 2018

1
Institute of Geography and Geology, University of Wuerzburg, D-97074 Wuerzburg, Germany
2
Alfred Wegener Institute, Helmholtz Center for Polar and Marine Research, Bussestr. 24, 27570 Bremerhaven, Germany
*
Author to whom correspondence should be addressed.
Remote Sens. 2019, 11(19), 2329; https://doi.org/10.3390/rs11192329
Submission received: 25 July 2019 / Revised: 26 September 2019 / Accepted: 6 October 2019 / Published: 8 October 2019
(This article belongs to the Special Issue Remote Sensing of Permafrost Environment Dynamics)

Abstract

:
Air temperatures in the Arctic have increased substantially over the last decades, which has extensively altered the properties of the land surface. Capturing the state and dynamics of Land Surface Temperatures (LSTs) at high spatial detail is of high interest as LST is dependent on a variety of surficial properties and characterizes the land–atmosphere exchange of energy. Accordingly, this study analyses the influence of different physical surface properties on the long-term mean of the summer LST in the Arctic Mackenzie Delta Region (MDR) using Landsat 30 m-resolution imagery between 1985 and 2018 by taking advantage of the cloud computing capabilities of the Google Earth Engine. Multispectral indices, including the Normalized Difference Vegetation Index (NDVI), Normalized Difference Water Index (NDWI) and Tasseled Cap greenness (TCG), brightness (TCB), and wetness (TCW) as well as topographic features derived from the TanDEM-X digital elevation model are used in correlation and multiple linear regression analyses to reveal their influence on the LST. Furthermore, surface alteration trends of the LST, NDVI, and NDWI are revealed using the Theil-Sen (T-S) regression method. The results indicate that the mean summer LST appears to be mostly influenced by the topographic exposition as well as the prevalent moisture regime where higher evapotranspiration rates increase the latent heat flux and cause a cooling of the surface, as the variance is best explained by the TCW and northness of the terrain. However, fairly diverse model outcomes for different regions of the MDR (R2 from 0.31 to 0.74 and RMSE from 0.51 °C to 1.73 °C) highlight the heterogeneity of the landscape in terms of influential factors and suggests accounting for a broad spectrum of different factors when modeling mean LSTs. The T-S analysis revealed large-scale wetting and greening trends with a mean decadal increase of the NDVI/NDWI of approximately +0.03 between 1985 and 2018, which was mostly accompanied by a cooling of the land surface given the inverse relationship between mean LSTs and vegetation and moisture conditions. Disturbance through wildfires intensifies the surface alterations locally and lead to significantly cooler LSTs in the long-term compared to the undisturbed surroundings.

Graphical Abstract

1. Introduction

Arctic landscapes have experienced rapidly increasing air temperatures of 0.6 °C per decade over the last 30 years, which is in an order of magnitude twice as high as the global average [1]. In particular, Arctic river deltas are considered to be majorly affected by rising temperatures, as they are located at the interface of the marine and terrestrial ecosystems [2]. These regions are extensively underlain by permafrost and therefore are sensitive to alterations in the thermal regime. As a result, large-scale environmental changes have occurred, including active-layer thickening, the formation of thaw slumps, progressive coastal erosion, and thermokarst development leading to surface subsidence, often followed by lake drainage [3,4,5,6]. Furthermore, the permafrost environments of the high latitudes store large amounts of soil organic carbon, which have the potential to further accelerate global warming after being released into the atmosphere as methane and carbon dioxide [7]. Yet the rising temperatures not only affect the state of the permafrost but also contribute to hydrological alterations, as well as vast land cover and vegetation changes [8]. Accordingly, a trend towards an overall greening of the Arctic has been observed where shrubs and woods are expanding northwards at the expense of tundra vegetation. This, in turn, leads to a decrease in surface albedo and, therefore, to increased uptake of energy by the land surface [9,10,11,12]. The environmental change of the Arctic occurs on multiple scales and in remote regions, hence remote sensing provides a unique opportunity to detect and to monitor these landscape dynamics in a continuous manner. Satellite-derived Land Surface Temperature (LST) characterizes the land–atmosphere exchange of energy and depends on a variety of surficial properties, such as vegetation type, soil, and plant moisture, or surface roughness [13]. LST may change with the alteration of the surficial properties, allowing environmental change to be characterized by means of time series remote sensing.
Previous studies have focused on the use of temporally dense but spatially coarse resolution remote sensing data (e.g., MODIS and AVHRR). AVHRR has been used to detect LST warming trends of approximately +0.72 °C per decade for the pan-Arctic regions between 1981 and 2005 [14,15]. MODIS derived LST has proven to be in accordance with in situ observed air temperatures and has been successfully applied to detect temperature trends, anomalies, and surficial changes [13,16,17,18]. Furthermore, by analyzing the relationship between the Normalized Difference Vegetation Index (NDVI) and LSTs on a pan-Arctic scale, the influence of the land surface thermal properties on arctic vegetation types and abundance has been revealed [19]. Nevertheless, Arctic landscapes are heterogeneous, and, therefore, a variety of processes cannot be detected using coarse resolution remote sensing data [2,16]. Additionally, the water-body fraction in the instantaneous field of view of the sensor can significantly lower the accuracy of the derived LST and represents an integrated signal of land and water surfaces [16]. The United States Geological Survey (USGS) Landsat Global Archive is able to compensate for these issues as it comprises several decades of optical data at relatively high spatial resolution.
Accordingly, this study presents the analyses of Landsat-derived LSTs and multispectral indices using the cloud computing capabilities of the Google Earth Engine (GEE). The time series analysis (1985 to 2018) was carried out for the Arctic Mackenzie Delta Region (MDR) in Canada. The MDR has been subject to major environmental changes with strongly increasing air and ground temperatures, and previous studies have revealed a large-scale greening of the landscape mostly expressed in shrub proliferation [20,21,22,23,24]. Exemplarily, the mean annual ground temperatures in the region have increased by 1–3 °C since the mid-1960s, although the magnitude of the trend is dependent on site-specific conditions [23,25]. Previous studies have analyzed land cover changes and trends in the MDR [5,26] and Lena delta [2] by means of time series regression analysis based on Landsat multispectral indices. Nevertheless, we lack knowledge of the long-term conditions of the surface thermal regime in the MDR at high spatial detail and how it is influenced by the physical surface properties, including vegetation cover, moisture regime, and topography. In periglacial environments, heat fluxes between the atmosphere and the ground are influenced by the presence of frozen sediments, and, in summer, a significant amount of the available energy is bounded to the melting of the active layer [27]. The ability to capture the development of summer LST detailed in space and time has, therefore, high relevance, as the temporal development of the LST may serve as a proxy on the recent active layer/permafrost development in the MDR.
Therefore, in this study, Landsat-based LSTs are retrieved by implementing a Single-Channel algorithm in the GEE. Long-term means of the LST, Tasseled Cap (TC) components, NDVI, and Normalized Difference Water Index (NDWI), as well as morphometric terrain features derived from the TanDEM-X digital elevation model (DEM), were calculated to represent different physical surface properties, including vegetation state, moisture regime, and relief situation. Firstly, the aim of the study was to investigate the influence of the surficial properties on the long-term mean of the surface thermal regime. Accordingly, correlation and multiple linear regression analyses were applied to establish a statistical relationship between the mean LST and surface features. Secondly, general temporal dynamics, as well as selected local land surface alterations including wildfire disturbance, were analyzed using the Theil-Sen regression slopes of LST, NDVI, and NDWI, which indicate the direction and magnitude of change over time between 1985 and 2018.

2. Materials and Methods

2.1. Study Area

The study area is situated in the continuous permafrost environment of the Mackenzie Delta Region in northern Canada between 67° to 70° N and 132° to 138° W (Figure 1a). The Mackenzie Delta is the second largest Arctic river delta covering an area of roughly 13,000 km2 [25]. The river flows northwards and, while bound by the Richardson Mountains in the west and the Caribou Hills in the east, diverges into several meandering channels that empty into the Beaufort Sea [25,28]. The region is situated at the transition of the boreal forest and the subarctic tundra biome, gradually divided by the tree line [25]. Accordingly, the landscape is diverse in terms of vegetation abundance and species, permafrost distribution, and the presence of surface water. The region has been subject to major environmental changes mainly linked to increasing ground and air temperatures [20,21,22,23]. In this study, four subregions (Figure 1a) were selected for the analysis that exemplarily highlight the ensemble of different landforms and land cover types along a north–south stretched gradient from the coast to the mountains. All regions were visited during fieldwork campaigns in the years 2012, 2013, and 2018.
The western part of the study area is located at the border of Yukon and the Northwest Territories, where the Richardson Mountains rise up to more than 1700 m above sea level, and terrain ruggedness is higher than in the rest of the landscape (Subregion 4). There, land cover is dominated by low tundra plant formations and extensive patches of dwarf shrubs while tall shrubs are rather seldom but can be found in wind-sheltered positions, whereas exposed hilltops, top slopes, and shoulders are only sparsely vegetated and often only covered by lichens [25,29]. In the center of the study area lies the delta itself, characterized by flat terrain and numerous lakes and channels (Subregion 3) [25]. The delta can be grouped into three major ecological zones [29]; spruce forest communities that established on sites less influenced by annual flooding dominate the southernmost part. The second zone is the transition between the two biomes, characterized by the increasing domination of willows and alder. Lastly, shrubs and herbs (willows and sedges) populate the tundra landscapes in the north. At the estuary, extensive wetland complexes, sand bars and islands have formed [25,29]. Adjacent to the delta in the southeast, a mosaic of open spruce forests and peat plateaus in the uplands dominate the landscape [25,30]. They are accompanied by tall shrubs that further north towards the uplands of the Caribou Hills change into dwarf shrubs consisting of willows, alder, and birches decreasing in size with increasing latitude (Subregion 2) [30]. Well-drained areas are populated by grasses and shrubs, whilst sedges dominate at moister locations [25]. In the northeastern part of the study area are the lowlands of the Tuktoyaktuk Peninsula as well as Richards Island, an outlier separated from the mainland by the East Channel (Subregion 1). These regions are characterized by rather subtle topography, rolling hills, and numerous depressions with thermokarst lakes, which continuously undergo environmental change expressed in thaw slump formation, lake drainage, or lake expansion [6,21,31,32].
The climate of the Mackenzie Delta Region is characterized by its pronounced seasonality and climatic gradients, which are determined by latitude, elevation, and coastal proximity—in particular, the presence of sea ice [25]. The mean annual air temperatures are –8.2 °C at Inuvik and –10.1 °C at Tuktoyaktuk for the period from 1981 to 2010. However, strong positive air temperature trends have been observed throughout the entire area, which is in accordance with the generally larger trend magnitude for high latitude regions and the Western Arctic of North America in particular [25,28,33]. Generally, warming trends seem to be strongest in autumn and winter and lowest during spring and summer, which is in accordance with pan-Arctic observations [25,33]. As a consequence, the mean annual ground temperatures in the region have increased by 1–3 °C since the mid-1960s [23,25]. The warming is believed to have caused widespread greening of the Arctic tundra landscapes expressed in large-scale shrub proliferation [24,26]. This is accompanied by an albedo reduction of the surface, amplifying further warming of the near-surface ground [24]. Additionally, the frequency of wildfires and the area affected by wildfires have increased, which has also been attributed to the observed temperature rise [34,35,36]. The potential consequences include permafrost degradation causing thermokarst development and active-layer thickening as well as vegetation alteration expressed by a distinct expansion of shrubs [22,24,37]. Whilst the overall temperature trend has caused widespread greening, wildfires can accelerate shrub expansion rather locally by more than double compared to unburned areas [24].

2.2. Data

The data used and the processing applied in this study are based on the cloud computing capabilities of the Google Earth Engine (GEE), which provides the opportunity for large-scale geospatial analysis [38]. GEE offers access to a variety of freely available archives of remote sensing data, among them the U.S. Geological Survey (USGS) Landsat Global Archive. Imagery is provided as raw digital numbers (DN) representing scaled radiance, calibrated Top–of–Atmosphere (TOA) reflectance, as well as surface reflectance (SR) data. All available Landsat-5, Landsat-7, and Landsat-8 images of all three processing types acquired from July to August between 1985 and 2018 with a maximum land cloud cover of 60% were included in this study, resulting in a total of 1699 scenes (Figure 2b). The months of July and August represent the peak growing season and have also been used in previous studies of Arctic landscapes, which allows for better comparability [2,9,26]. The high latitude of the study area results in a strong overlap of the WRS-2 paths, which increases the acquisition frequency and thus resulted in an overall dense time series at the pixel level (Figure 2a).
The Landsat data is complemented by topographic parameters derived from the high-resolution digital elevation model (DEM) of the TanDEM-X mission provided by the German Aerospace Center (DLR) (see Acknowledgments) [39,40,41,42]. The DEM covered most of the region of interest (Figure 1c) and had a pixel spacing of about 12 m, after reprojecting to UTM Zone 8N using the WGS1984 ellipsoid. The land surface parameters slope and aspect were processed using the DEM, and the aspect was then transformed using a sine-function to avoid circular data and retrieve the STAGE parameter—the northness of the terrain exposition. Additionally, the potential solar insolation in kWh/m² was estimated in SAGA GIS (http://www.saga-gis.org/) according to the approach of [43,44]. Further, the DEM was used in hydrographic modeling; in the pre-processing, a highly detailed governmental vector dataset on the hydrography (open.canada.ca), including information on ocean, lakes, ponds, rivers, and channels, was applied. Using the pre-processed DEM, the catchment size of each pixel was processed using a Multi-Flow-Direction approach (MFD) [45]. Finally, the Topographic Wetness Index (TWI) was calculated [45]. This index displays the logarithmic ratio between the size of the catchment and the local slope: higher index values characterize flat regions with rather large catchments, whereas low index values indicate steep locations with rather small catchments.

2.3. Methods

Figure 3 illustrates the processing chain that was applied in this study, which includes processing the Landsat data and retrieving the LST as well as the multispectral indices, in order to calculate the statistical temporal metrics. The GEE implementation for the retrieval of the temporal statistical metrics, including the main feature of Landsat LST derivation, is available on GitHub (https://github.com/leonsnill/lst_landsat).

2.3.1. Pre-Processing and Retrieval of Multispectral Indices

The Landsat imagery was masked for clouds, cloud shadows, and snow or ice using the Quality Assessment Band of the SR product generated from the CFMASK algorithm. The SR product was also used to obtain the NDVI (Equation (1)) and NDWI (Equation (2)) for each image [46,47]:
N D V I = ρ N I R ρ R E D ρ N I R + ρ R E D   ,
N D W I = ρ N I R ρ S W I R 1 ρ N I R + ρ S W I R 1   ,
where ρ x is the reflectance in the corresponding part of the electromagnetic spectrum. The NDWI based on [47] and used here focuses on the water content in vegetation rather than water bodies as the same-titled index by [48]. NDVI and NDWI were chosen as they are sensitive to chlorophyll content, vegetation water content, as well as subpixel water-fraction and subpixel vegetation-fraction, respectively [49]. However, these indices tend to saturate in high-density canopies and reduce available variance by disregarding parts of the spectral feature space. Contrary, the TC components preserve variance in the data while allowing for reducing the overall dimensionality of the data and have been used to detect environmental changes in the Arctic [2,26,50]. Accordingly, the TOA product was used to derive the Landsat-specific Tasseled Cap transformations greenness (TCG), brightness (TCB), and wetness (TCW) (Equation (3)):
T C x =   ρ B C B + ρ G C G + ρ R C R + ρ N I R C N I R + ρ S W I R 1 C S W I R 1 + ρ S W I R 2 C S W I R 2 ,
where the sensor and band-specific coefficients Cx that were used in this study are summarized in [2], and ρ x are the reflectance values in the corresponding parts of the electromagnetic spectrum.

2.3.2. Retrieval of Land Surface Temperature

For Landsat thermal infrared data with a certain channel width, one can obtain an effective at-satellite temperature BTsen (Equation (4)) based on an approximation of Planck’s law as follows [51]:
B T s e n = K 2 l n ( K 1 L s e n + 1 ) ,
where K1 and K2 are band-specific thermal conversion constants provided with the metadata and Lsen refers to the spectral radiance in W/(m2⋅sr⋅μm), which can be obtained by applying the band-specific rescaling factors Gain and Offset also provided with the metadata file to the pixel values (DN) (Equation (5)):
L s e n = G a i n D N + O f f s e t
The SR collection directly provides the calibrated brightness temperature (BT) needed for the retrieval of the LST, as described in the following section. This also includes at-sensor radiance (L) of the corresponding thermal band, which was calculated by applying the radiance rescaling factors provided in the metadata file to the DNs of the raw L1TP collection [52].
In this study, a Single-Channel (SC) algorithm was used to retrieve LSTs that requires knowledge of the surface emissivity and the state of the atmosphere. This method was chosen as it is comparably simple to implement if the parameters are known, and because it has proven to be accurate and applicable for sensors with only one thermal band, such as in the case of Landsat TM and ETM+ [53,54]. Landsat-8 offers two bands in the thermal infrared region allowing for the use of a Split-Window (SW) algorithm. These are widely used and have proven to be more accurate than SC methods, but TIRS has been subject to contamination by stray light, especially in band 11 [55]. Therefore, it is advised not to implement an SW algorithm as it may lead to higher uncertainties in the retrieval of LST [56]. For that reason and for the sake of comparability with TM and ETM+, the generalized SC algorithm developed by [54] was used to derive LST for all three sensors. Accordingly, on the basis of Planck’s law and a radiative transfer model, the LST (Equation (6)) can be retrieved as follows [54,57]:
L S T = γ [ 1 ε ( ψ 1 L s e n   +   ψ 2 )   +   ψ 3 ] +   δ ,
where Lsen is the spectral radiance at the sensor in W/m2⋅sr⋅μm (Equation (5)), ε is the surface emissivity, ψ are so-called atmospheric functions (AFs), and γ (Equation (7)), as well as δ (Equation (8)), are parameters based on Planck’s law:
γ = B T s e n 2 b γ L s e n ,
δ = B T s e n B T s e n 2 b γ ,
where BTsen is the brightness temperature from (Equation (4)), and bγ is a sensor-specific constant taking a value of 1256 K, 1277 K, or 1324 K in the case of TM, ETM+, or TIRS, respectively [57,58]. The AFs describe the state of the atmosphere with regards to transmissivity, upwelling, and downwelling radiation and are approximated versus the atmospheric water vapor (WV) content using a second-degree polynomial fit (see [22] for details):
[ ψ 1 ψ 2 ψ 3 ] = [ c 11 c 12 c 13 c 21 c 22 c 23 c 31 c 32 c 33 ] [ w v 2 w v 1 ]   .
The coefficients cij are retrieved by simulation using different atmospheric soundings databases, resulting in different coefficients for each sensor. The coefficients used in this study can be found in Appendix A (Table A1) and are best suited for high latitude environments with usually low WV content [57].
The atmospheric WV content was retrieved within the GEE based on reanalysis data from the National Centers for Environmental Prediction (NCEP) and the National Center for Atmospheric Research (NCAR) [59]. The NCEP/NCAR Reanalysis Project provides the total column water vapor at a global scale with a spatial resolution of 2.5 arc degrees and with a 6-hourly analysis field (00:00, 06:00, 12:00, 18:00 UTC). Accordingly, the WV content estimation closest to the observation time of each scene was selected rather than taking the mean of the two closest estimations, as the vast majority of Landsat scenes were close to the actual 18:00 UTC analysis field. Reanalysis data is considered to be accurate for the atmospheric correction of the thermal infrared data, and NCEP/NCAR WV data, specifically, has proven to yield accurate results when retrieving LST from Landsat data [60,61].
The land surface emissivity was derived for each time step using the Simplified NDVI Threshold Method (SNDVITHM) described by [62]. Applying this method, emissivity is obtained from the NDVI, as both parameters show a linear relationship, and the emissivity for soil and vegetation is almost consistent within the 10 μm to 12 μm range of the electromagnetic spectrum [60,62,63]. Therefore, when the study area mainly consists of bare ground and vegetation cover, this method has proven to be an accurate estimator [63]. For this method to be implemented, one must choose certain threshold values representing the emissivity and index values of full vegetation cover (εv and NDVIv), as well as bare soil (εs and NDVIs). The pixel emissivity values may then be derived accordingly [62] (Equation (10)):
ε = { ε s   ε s + ( ε v     ε s ) P V           ε v N D V I < N D V I s N D V I s   N D V I     N D V I v N D V I >   N D V I v
with PV being the fraction of vegetation cover [64] (Equation (11)):
P V = ( N D V I N D V I s N D V I v N D V I s ) 2 .
NDVIs and NDVIv were chosen to be 0.2 and 0.6, respectively. The emissivity values εs and εv were set to 0.97 and 0.985, respectively. The latter were selected based on calculations from previous studies and emissivity curves in the thermal infrared region [60,63,65,66,67]. Furthermore, to account for waterbodies, NDVI values below zero were assigned an emissivity value of 0.99. Finally, LST was calculated for each scene using (Equation (6)).

2.3.3. Statistical Analysis

The final processing included mosaicking the multi-temporal layer stacks and correcting them for outliers, which were found to be mainly introduced by insufficient cloud and cloud shadow detection of the CFMask algorithm in certain scenes. The correction was achieved by only including pixels above a temperature threshold of 8 °C. Mosaicking scenes of the same date, hence the Landsat WRS-2 path, became necessary as the overlap between adjacent WRS-2 rows introduced errors when calculating the standard deviation, as well as the T-S slope parameter. This could not be attributed to a particular issue, yet we believe it is related to the increased number of observations in the overlapping part of the scenes. After that, the statistical temporal metrics for all parameters (LST, NDVI, NDWI, and TC), namely the mean, the standard deviation, and the slope coefficient of the Theil-Sen regression, were calculated. A tabular summary of all investigated parameters is provided in Appendix B (Table A2).
In order to quantify the relationship between the surface thermal regime and certain surface properties (e.g., vegetation cover, surface moisture, topography), correlation and regression analyses were conducted using the long term means of the multispectral indices and topographic features as explanatory variables of the mean LST. This was done for the entire study area, as well as the chosen subregions depicted in Figure 1. The statistical analysis was carried out in Python using the statsmodels module. In a first step, Pearson’s linear correlation coefficient (R) was used to detect statistical relationships between variables, which was followed by the interpretation of the causality of these relationships. Secondly, Pearson’s R served as a proxy in the initial selection of variables used in the multiple linear regression models to only include correlated variables to the mean LST (here R > 0.3). As the correlation matrices revealed correlations among independent variables, the Variance Inflation Factor (VIF) was calculated for each variable in each set of selected variables of each region (i.e., study area and four subregions) to account for multicollinearity. In this study, the widely used threshold of VIF < 5 was used to select the variables. The selected variables for each model region were then used to build the regression models to assess the relative importance of each factor on the mean LST and to assess if multi-variable models increase the predictive ability.
Secondly, the slope coefficient of the T-S regression was used to detect the direction and magnitude of change over time to reveal general trends and dynamics of the land surface. Accordingly, the pixel-based slopes of LST, NDVI, and NDWI were analyzed in their overall spatial expression as well as their spectral-temporal behavior for selected plots in the subregions. The overall pixel-based slopes were calculated within the GEE using the sensSlope-Reducer function whilst the selected single plots were calculated in Python using the theilslopes function of the scipy module. The T-S regression has proven to be a robust estimator of trend direction and magnitude, being insensitive to up to 30% of outliers, and has already been successfully applied to detect spatio-temporal landscape dynamics in an Arctic environment [2,5,26,50,68,69]. It is calculated as the median of all slopes between every pair of given values [70].

3. Results

3.1. Mean Summer Land Surface Temperature, NDVI, and NDWI

Figure 4 shows the processing results of the spectral-temporal metrics mean and standard deviation of the summer LST, NDVI, and NDWI for the entire study area. Waterbodies were masked using the highly detailed governmental GIS-database on the hydrography of Canada. For the entire study area, most of the observations (>98%) of mean summer LSTs range between 15 °C and 28 °C However, distinct spatial patterns arise, and the distribution was bimodal with a smaller peak at approximately 18 °C and a second larger peak at approx. 22.5 °C. The first peak corresponds to the cooler Mackenzie Delta, which exhibits the lowest LSTs with values mostly less than 20 °C. The outer delta (towards the western estuary and the Shallow Bay), the north of Richards Island, as well as the northern coastal area of the Tuktoyaktuk Peninsula show moderate mean summer LSTs with the majority of values ranging between 20 °C and 22 °C. The uplands of the Caribou Hills, of the Anderson Plain, the Peel Plateau, the Richardson Mountains, as well as the slopes bounding the Mackenzie Delta in the west, are characterized by mean summer LST values above 22 °C. However, incised valleys and north-facing slopes of the Richards and the Richardson Mountains show LST values frequently below 20 °C. Noticeably, for the Anderson Plain and the Caribou Hills, clear and distinct patches of lower and higher LSTs are observed. These locations correspond to the extent of former wildfires; hence, the disturbances due to fire events are expressed in the mean LST. Most of these regions are further characterized by a higher standard deviation of the LST (Figure 4d). The distribution of the standard deviation of the summer LST showed unimodal symmetric frequency distribution, and deviations are of 4.5 °C magnitude on average. The majority of the deviations (98%) range between 3 °C and 6.5 °C.
Figure 4b displays the processing results of the mean NDVI, clearly highlighting regions of sparse and dense vegetation cover. Exemplarily, the uplands of the Richardson Mountains and the outer Mackenzie Delta are characterized by lower NDVI values, whereas the highest NDVIs are observed for regions in sheltered positions, e.g., valleys and mid-slopes, but also for regions that were formerly affected by wildfires. The standard deviations of the summer NDVIs (Figure 4e) are rather small, and most of the observations (98%) display variations less than 0.1. Most invariant are regions that seem not to be covered by vegetation. The distribution of the mean NDVI showed that all values are positive and 98% of all observations are in the range between +0.2 and +0.8. The frequency distribution is centered at a value of approximately +0.65. Overall, deviations are small, and 98% of the standard deviations are less than 0.012 and range between 0.005 and 0.008.
The results of the processing of the NDWI are displayed in Figure 4c. Negative values of the mean summer NDWI are observed over the outer Mackenzie Delta and the Richardson Mountains. NDWI values closely around zero are found for the Anderson Plain, northern Richards Island, and the northern coastal areas of the Tuktoyaktuk Peninsula. In contrast, positive NDWI values are most frequent and are observed for most areas of the Mackenzie Delta, the inland of the Tuktoyaktuk Peninsula, and especially for some of the regions formerly affected by wildfires. The majority of the values (98%) range between +0.5 and −0.5. The variations of the NDWI in time are small and 98% of all observations are characterized by standard deviations less than 0.006 (Figure 4e). The strongest variations in time are found for the Anderson Plain and the outer Mackenzie Delta. The histograms of the NDWI features show that 98% of the region displays index values between −0.4 and +0.4. The frequency distribution is strongly left-skewed, and the majority of the values are greater than zero. The distribution is centered at an NDWI index value of +0.18. As mentioned above, variations over time are of low magnitude, and 98% of all observations are in the range from 0.0025 to 0.012. Overall, the patterns of mean LST, NDVI, and NDWI match well where drier and less densely vegetated areas generally exhibit higher LSTs and vice versa.
The LST at the site “East Channel” clearly shows differences in the mean summer LST between the incised channel and the uplands of Richards Island and the Tuktoyaktuk Peninsula. LST values of the lower elevated terrain with closer proximity to the channel are, on average, approximately 1–2 °C cooler than the LST values of the uplands. The top slopes of the cliffs show higher LST and reduced NDVI and NDWI values. Field visits proved that these locations were more exposed and exhibited lower vegetation coverage than the adjacent sheltered leeward sites. The results of the study area “Inuvik” indicate the influence of two former fires on the mean summer LST, NDVI, and NDWI. The largest fire event occurred in 1968, west of Lake Noel [71,72]. The former extent of the fire is clearly visible in the data, and the mean summer LST is 1–2 °C colder than the undisturbed surroundings. Similarly, the fire event of 1968 is clearly visible in the mean of NDVI and NDWI, as the area is characterized by higher index values indicating a densely vegetated and moister environment. Similarly, the relation of the features LST, NDVI, and NDWI is deduced from the study area “Delta”. The highest LST values are found in the west of the area, off the main channel. Here, the higher values in LST are associated with moderately high NDVI and NDWI values, yet those are the lowest within this subplot. In contrast, the site “Richardson Mountains” showed a clear relation between the LST and terrain exposition (STAGE) as northward facing slopes are significantly colder than southward facing slopes with differences of up to 8 °C.

3.2. Correlation and Regression between LST and Environmental Factors

The relationship between the features was quantified in linear correlation and multiple regression analyses. Figure 5 displays a subset of the correlation matrices for each region showing Pearson’s R. The predicted-observed density plots of all regression models are shown in Figure 6. As described in Section 2.3.3, the features have been selected based on the correlation results and filtered by the VIF to account for multicollinearity as well as an individual assessment based on the model fit and complexity (number of independent variables) by taking into account the Bayesian information criterion (BIC) in the model selection process. Detailed information on the output results of the models are provided in Appendix C (Table A3, Table A4, Table A5, Table A6 and Table A7).
For the entire study area, the highest correlation is observed between the mean summer LST and the mean TCW (−0.48), as well as the northness of the terrain (STAGE) (−0.42). Furthermore, only the means of the NDWI and TCB also showed a weak to moderate correlation with the mean LST with coefficients of −0.32 and 0.32, respectively. The model using the mean TCW and STAGE as predictors explains 31% of the response variables variance with an RMSE of 1.73 °C (Figure 6a). In a test, all variables were used in the VIF selection process regardless of their correlation to the mean LST, and the resulting model consisted of eight explanatory variables (mean of NDVI, TCB, TCW, as well as DEM, TWI, FlowAcc, STAGE, and WaterDist), which only increased the explained variance by 4%. Additionally, by comparing the relative variable importance based on the z-score standardized coefficients, the mean TCW and STAGE can clearly be regarded as most influential on the mean LST and provide a comparably similar model fit using fewer parameters.
For the subregion “East Channel”, the highest correlations of the mean LST are observed with the mean TCW (−0.79) and TCB (0.72), as well as the DEM (0.51). Further, the mean NDWI (−0.46), the TWI (−0.41), the mean TCG (0.37), and the WaterDist (0.33) correlate with the target variable. Initially, the VIF analysis revealed the highest score for the mean TCW; however, given the higher correlation with the LST, the NDWI having the second largest score was removed instead. The correlation between TCW and NDWI (0.72) also indicated large parts of the variance of the NDWI being captured by the TCW. The final model shown in Figure 6b shows a great fit with an R2 of 0.73 and an RMSE of 0.81 °C. By comparing the standardized coefficient estimates (Table A3, Table A4, Table A5, Table A6 and Table A7), the mean TCW, the TWI, and WaterDist features revealed to be of the highest relative influence.
For the site “Inuvik”, terrain features are far less influential than the means of the multispectral indices that correlate best with the mean LST. Again, the mean TCW shows the strongest correlation with an R of −0.75, followed by the mean NDWI (−0.73). This is also the only subregion in which the mean NDVI exhibited a stronger correlation (−0.54) to the LST. Furthermore, the TCG and TCB were included in the selection process as they showed coefficients greater than 0.3. The VIF analysis yielded the model depicted in Figure 6c that is constituted of three parameters representing three different surface characteristics, namely the moisture/water (TCW), vegetation (NDVI), and soil properties (TCB). The model explains 63% of the variance with an RMSE of 0.83 °C. Again, the mean TCW shows the greatest influence on the target variable, followed by the mean NDVI.
The mean LST of the ”Delta” subregion is exclusively correlated with moisture/water features, namely the mean NDWI (−0.84), TCW (−0.69), and the TWI (0.42). The VIF did not indicate collinearity problems, and the model depicted in Figure 6d, therefore, consists of all three features. Almost three quarters of the variance (74%) of the mean LST can be explained with an RMSE of 0.51 °C, yet again, one parameter alone contributes to the largest part in explained variance, i.e., the mean NDWI.
In contrast to the previous areas, the mean LST in the “Richardson Mountains” subregion is strongly correlated with the topographic features STAGE (−0.66) and Insolation (0.61) and to a lesser extent with the means of TCB (0.54) and TCW (−0.40). Due to strong collinearity, the VIF suggested to exclude the Insolation parameter, and the model shown in Figure 6e is largely driven by the northness of the terrain (STAGE) and the mean TCB, whilst the TCW only shows little relative influence on the mean LST. The model explains approximately 48% of the target variables variance, and the RMSE of 1.61 °C reveals large deviations between the predicted the observed mean LST.

3.3. Temporal Dynamics of LST, NDVI, and NDWI

Between 1985 and 2018, the MDR is characterized by surface dynamics that vary across space in their direction and magnitude, as can be derived from the T-S trends of the LST, NDVI, and NDWI (Figure 7). A large-scale wetting and greening trend can be observed according to the NDVI and NDWI, respectively. Generally, this trend increases with latitude and is strongest in the coastal lowlands of Richards Island and the Tuktoyaktuk Peninsula, exhibiting large slope coefficients of 0.07 per decade and above. An additional gradient following altitude can be observed with decreasing trends of the NDVI and NDWI as elevation increases in the Richardson Mountains. For the entire study area, a mean decadal increase of 0.031 ± 0.020 (one standard deviation) and 0.027 ± 0.022 is observed for the NDVI and NDWI, respectively. Locally, patches that range from a few hundred meters up to dozens of kilometers either indicate strong greening and wetting or browning and drying, which mostly reflects areas affected by former wildfires (Figure 7, fire extents are highlighted in white). The alluvial plain in the central study area shows strong local variation in the direction of the slope coefficient, yet of mostly smaller magnitude compared to the surrounding lowlands. The large-scale exception is the browning and drying of the outer delta that is unique to the entire study area. In general, the T-S trends of the NDVI and NDWI revealed a moderate to strong correlation of 0.57, hence wetting and greening or browning and drying may often be accompanied by each other and vice versa.
T-S trend patterns of the LST show similarities to the T-S trends of NDVI and NDWI. Positive trends of the indices are mostly associated with a cooling of the land surface, whilst strong browning and drying, as present in the estuaries of the outer delta, is associated with an increase in surface temperatures. Accordingly, the correlations between the T-S trends of the LST with the NDVI (−0.23) and the NDWI (−0.49) are weak to moderately strong. Overall, LSTs exhibit a mean decadal decrease of −0.345 ± 0.527 °C in the MDR. Accordingly, the T-S trends are spatially diverse, strongly regionalized, and not unidirectional. The northern section of the study area, and especially the coastal highlands, exhibit strong cooling trends of approximately −0.5 °C to more than −1.5 °C per decade. In the lowlands, where numerous thermokarst lakes dominate the landscape, these cooling trends are less pronounced and, locally, even shift to positive trend slopes. The latter is especially pronounced in the northernmost coastal areas. The outer delta is again an exception, as positive T-S trends of the LST are associated with negative trends of the T-S trends of NDVI/NDWI. The alluvial plain is characterized by heterogeneous patches of positive and negative slope coefficients, while the latter are more frequent and of greater magnitude. The southwestern regions of the study area (towards the Richardson Mountains and the Peel Plateau) are dominated by positive T-S trends of 0.5 °C per decade on average, yet local extremes exceed 2 °C per decade. In contrast, the southeastern regions of the study area (towards the Anderson Plain) are characterized by mostly decreasing temperatures, in particular, areas associated with wildfires exhibit a strong cooling trend, whereas strong positive T-S trends of the LST are rather small-scale.

3.4. Local Temporal Dynamics of LST, NDVI, and NDWI

In addition to the analysis of large-scale surficial dynamics in the MDR, small-scale changes were investigated for two subregions with a size of 30 × 30 km. The first was the coastal “East Channel” region (Figure 8) that is characterized by fluvial processes, subtle topography with rolling hills, and depressions containing numerous thermokarst lakes. The region exhibits T-S slope coefficients in both directions (i.e., cooling and warming), whereas most of the land surface shows a cooling trend, which is particularly pronounced in the lake-rich lowlands in the southeastern part. Moderate to strong warming trends, on the other hand, are almost exclusively found in the fluvial environment of the East Channel and the low-lying regions of southern Richards Island. The NDVI and NDWI slope coefficients reveal an extensive wetting and greening of the land surface, with only a few exceptions in the surrounding of thermokarst lakes on Richards Island, as well as the fluvial islands of the East Channel. The lowlands associated with the strong LST cooling trend are, accordingly, those characterized by the most extensive and pronounced wetting and greening: at plot location e (Figure 8e), the NDVI and NDWI increased by more than 0.065 per decade between 1985 and 2018, with only little variance of the index values over time. The significant increase of the indices is also illustrated by the narrow upper and lower limit of the 95% confidence interval of the slope coefficient. The LST trend line indicates an associated strong cooling of the land surface; however, given the larger variance over time, the confidence interval of the slope coefficient is significantly larger.
The second location (Figure 8f) is situated on a fluvial island of the East Channel and represents an area of highly active hydrological dynamics, including flooding and the erosion and accumulation of sandbars. The temporal trajectories of NDVI and NDWI demonstrate these dynamics given their cyclic shape over the observation period with three peaks around the years 1988, 1999, and 2013. Although linear trends do not capture this cyclic behavior, they might indicate long-term developments towards a drier or wetter environment. Generally, the slopes of NDVI and NDWI indicate the development towards a greener and relatively drier environment, whilst the LST seem to be increasing. Plot g is located at a formerly drained lake characterized by an increase in vegetation cover and surface moisture. The trajectories of the indices are clearly increasing in a linear fashion with some variance and abrupt changes between individual years. The T-S slope indicates that LSTs have decreased substantially analogously to plot location e.
The second subregion, “Inuvik” is depicted in Figure 9, which illustrates the temporal dynamics analogously to Figure 8. This region has been subject to wildfires with two extents (1968 and 2012) being highlighted in Figure 9a. The undisturbed sites in the north are predominantly characterized by a steady wetting and greening of the landscape (plot e). The T-S trend of the LST indicates a subtle cooling of the land surface, accompanied by an increase in vegetation cover and moisture. It is important to note that Landsat-7 SLC-off patterns are visible in the subregion and create an overall heterogeneous and patchy LST image. The region of the 1968 wildfire depicted in plot location f was subject to an extensive cooling of the land surface, accompanied by a strong increase of the NDVI and NDWI of up to 0.3 between 1985 and 2018. The area affected by the most recent fire in 2012, in contrast, exhibits browning, drying, and warming trends. The disturbance is clearly visible in the temporal trajectories of the parameters as a sharp decline of the multispectral indices can be observed (Figure 9g, dashed line). Overall, the borders of the wildfire events are present in all three slope images.

4. Discussion

4.1. Processing of LST Using Dense Landsat Time Series

This study utilized a single channel algorithm for the processing of the LST from Landsat datasets in the GEE that relied on the framework of [54,58,66] and included products on the column water vapor from the NCEP/NCAR Reanalysis Project. Like other studies that have dealt with Landsat data for long investigation periods and large study areas, no absolute referencing of the LST products was feasible, as no ground-truth data existed. However, the LST retrieval approach from Landsat data is well established and can be considered sufficiently accurate [58,61]. The presented results, therefore, represent an authentic and plausible remote estimation of the LST for the MDR, following the recent state of practice that also includes automatic cloud masking and outlier removal.
Nevertheless, it should be noted that LST exhibits strong diurnal variations and is highly sensitive to short-term synoptic variations of air temperature and insolation. The variance of the LST is, thus, inherently higher than, for instance, information on vegetation cover at peak growing season, especially in Arctic environments. For the LST, this results in rapidly changing insolation rates and heating of the land surface. Furthermore, the presence of near-surface permafrost introduces additional uncertainties as the ground properties (like active layer thickness, air and ice content) are usually unknown and highly heterogeneous in space and time. Consequently, these factors influence the heat fluxes and the resulting LST. Furthermore, investigating Landsat-based LST development and associated heat fluxes over all seasons is not feasible and restricted to the summer. There is, therefore, a very short time window within which to capture Landsat acquisitions that are suited for the analyses of the LST, even though WRS-path overlaps increase the per-pixel data density (Figure 2).
However, calculating long-term means of the LST compensates for these issues and creates spatially extensive information on the mean summer LST. This is reasoned by the fact that the surface thermal regime is, on average, determined by the physical surface characteristics and not by atmospheric conditions [13]. Accordingly, this allows for characterizing the surface thermal regime of the Arctic MDR by revealing patterns of cold- and hotspots and influencing factors that may serve as a proxy for permafrost distribution and development. Furthermore, the high spatial detail at the pixel level of LSTs derived from Landsat (resampled 30m) considerably improves the analysis in heterogenous permafrost environments, whilst using MODIS LST data (approximately 1000 m) may limit the assessment of local conditions [13].

4.2. Relation of LST to Other Environmental Variables in the Mackenzie Delta Region

The correlation and regression analysis allows for establishing a statistical relationship-basis between the mean LST and the surface properties expressed by the TC components, the NDVI, and NDWI, as well as the topographic features derived from the TanDEM-X DEM. Former studies have identified the NDVI and NDWI to be sensitive to the content of chlorophyll (NDVI), the moisture regime, and the vegetation water content, respectively (NDWI) [47,49,73]. Although causality cannot be derived directly from the pure correlation between two parameters, the analysis revealed convincing dependencies between the thermal regime, the vegetation cover, moisture regime, and topographic situation.
From a model perspective, the plant and soil moisture regimes of the landscape represented by the mean TCW generally explain the largest proportions of the variance in the mean LST. Increasing moisture leads to decreasing surface temperatures, which can be attributed to increased evapotranspiration rates and associated latent heat fluxes [13,27]. The delta itself with its numerous channels and the vast wetland complexes of the entire study area provide great moisture supply and are therefore highly prominent in the data with corresponding importance in the regression models. These findings are in accordance with [13] who found wetlands to exhibit the coolest temperatures among different land cover types in an Arctic environment.
However, the models of the different subregions were substantially better in their model fit and predictive ability compared to the entire study area. The diverse environmental setting of the MDR with mountainous regions in the west and south, the delta, and lowlands, as well as the extensive flat northern and eastern coastal areas, resulted in overall low correlation coefficients of the features. Accordingly, no single variable could explain the spatial patterns of the mean LST directly for the entire area on a high level of determination. Overall, the second most important variable was the northness of the terrain (STAGE), which can be explained by the vast stretching Richardson Mountains and rolling hills where southern facing slopes receive significantly more radiation. Furthermore, these areas are inherently well-drained, giving rise to generally drier surface conditions with increased sensible heat fluxes. On the contrary, northern facing slopes—although possibly equally well-drained—exhibit much colder LSTs, as can be seen in the long-term thermal mean (Figure 4a). As these areas constitute the majority of extreme cold- and hotspots in the MDR, the resulting feature space is well suited for predicting mean LSTs in the regression models, which is particularly true for the homogenous subregion “Richardson Mountains”.
On the contrary, for the subregions “East Channel”, “Inuvik”, and “Delta”, it was found that morphometric features offered only a little information on the spatial variability of the mean LST. This is probably because the relief is rather subtle and differences in the exposition are less pronounced. Here, the influence of the different vegetation cover and associated moisture regimes was more important in explaining spatial LST variations. In comparison to the moisture indices, the vegetation indices NDVI and TCG generally correlated less with the mean LST. This may be attributed to two factors: firstly, LSTs are largely controlled by evapotranspiration rates and they are better captured by moisture indices like the TCW. Secondly, although transpiration is controlled by the vegetation, different vegetation types (deciduous vs. coniferous) and the leaf area index are of greater importance [27]. However, this may not be captured by the indices, as the NDVI, for instance, tends to saturate in high-density canopies. Accordingly, information on the vegetation type and land cover would serve as a valuable explanatory variable when trying to understand the spatial expression of the mean LST [13].
Overall, large portions of the target variable’s variance remain unexplained. The regression results, therefore, suggest that the mean LST needs to be explained by multiple variables that capture a variety of surface and ideally, sub-surface characteristics. The reason for this might be the influence of different soil properties, including the thermal conductivity of the soil, soil moisture, active layer thickness, proximity to the permafrost table, and permafrost temperature [27]. These factors have a direct influence on the LST, but may also indirectly influence the LST, since vegetation cover and moisture regimes are dependent on sub-surface conditions. Furthermore, in permafrost environments, heat fluxes between the atmosphere and the ground are heavily influenced by the presence of frozen sediments: as illustrated by [27], in spring and summer, a significant amount of the available energy is bounded to the melting of the active layer and the permafrost, respectively, while taking a degradation into account.
The study focused solely on the land surface by not including (permanent) water bodies in the analysis, which is an important consideration as the large contrast in LSTs between water bodies and land surfaces can result in a pronounced distinction in the feature space: for instance, the mean NDVI may correlate highly with the mean LST when including water bodies due to the distinct clusters of water and non-water in the feature space, which may strongly limit the interpretability and understanding of the relationship between the variables and ultimately the potential to model LSTs.

4.3. Temporal Changes

The extensive greening of the land surface, as indicated by the T-S slope of the NDVI, is in accordance with previous studies in the MDR and on the pan-Arctic scale, which found that most tundra landscapes have been subject to rapid and vast greening [74]. Nitze & Grosse [2] observed the strongest vegetation trends in the Arctic Lena delta in coastal proximity, and Fraser et al. [26] found that the coastal areas of the Tuktoyaktuk Peninsula were most extensively affected by greening processes in the MDR between 1985 and 2011. This study reveals that these trends have continued to persist in the MDR for the additional timeframe observed and despite the aforementioned limitations of the LST trend product, have been accompanied by a cooling of the land surface, indicating the associated changes of the thermal regime with the alteration of the vegetation. Nitze & Grosse [2] attributed the strong increase in vegetation indices in coastal proximity to the rapid decline in sea ice cover over recent decades. In fact, the beginning of sea ice melt in the Beaufort Sea has exhibited a large negative trend of more than 10 days per decade, which is amongst the most rapid declines in the entire Arctic [75]. As sea ice concentration is a major controlling factor of the thermal regime, air temperatures have increased accordingly, causing an overall greening of the landscapes [26,71,74].
Local surficial changes over time can be well studied using time series of multispectral indices and LST. Generally, the local T-S trends in the East Channel and Inuvik subregions reveal the overall picture observed in the MDR, which is mostly characterized by greening and wetting processes. On the contrary, surface dynamics, including fluvial activity and wildfires, either reverse or enhance these trends, and this is also where the largest slope coefficients are present in the time series. Initially, fire events lead to a strong albedo reduction as the surface is charred and the soil organic layer degraded, which manifests in drier and warmer surface conditions with increased sensible heat flux [27]. A deeper heat flux into the ground causes near-surface ground temperatures to rise and active-layer thickness to increase [24,37]. This is followed by secondary succession, where higher soil moisture contents, due to the thaw and increased availability of nutrient matter, lead to an intense regrowth, and strong greening and wetting trends can be observed. Whilst the overall temperature trend has caused widespread greening, wildfires can accelerate shrub expansion rather locally by more than double compared to unburned areas [24]. Accordingly, the 1968 wildfire area that was formerly dominated by coniferous vegetation, mosses, and lichens became subject to fast and intensive shrubification by deciduous species and an expansion of grasses [75,76]. This systematic was observed in the data of the Landsat time series features comparing the 2012 and the 1968 fire events, which occurred in what were originally almost identical environmental settings. The 1968 wildfire area is characterized by lower mean LSTs than the undisturbed surroundings, while mean NDVI and NDWI values were usually higher and have increased strongly between 1985 and 2018. Landsat summer mean LST indicated that the 1968 region was on average 1 °C to 2 °C cooler than the undisturbed surroundings, which may be attributed to the higher soil moisture and correspondingly increased evaporative cooling, as well as to changes in vegetation species increasing evapotranspiration. This underpins the fact that moisture indices in the “Inuvik” subregion were highly correlated with the mean LST and further reinforces the necessity to include information on the land cover or vegetation types in the regression models. In contrast to the 1968 wildfire region, the area of the 2012 wildfire event exhibited elevated temperatures and sharp declines in the index values, which fits the description of Eugster et al. [27] and indicates an initially reduced albedo due to the increased soil signal. Increased wildfire disturbance, therefore, accelerates the process of greening and moisture supply leading to cooler LSTs.
The linkage between NDVI, NDWI, and LST is also expressed in the temporal domain of undisturbed sites (i.e. not affected by wildfires). The examples on the subsets “East Channel” and “Inuvik” highlighted the inverse relation of the variables where surface cooling was observed along with wetting and greening. These results indicated, as well, that the trend analysis of NDWI and NDVI can be performed on a higher level of confidence as temporal variations were of lower magnitude than the temporal variations of the LST. This may be because, from a processing perspective, the number of observations per pixel seemed to create patterns in the LST T-S trend product within the GEE. The scene boundaries are also visible in the data, and because of the aforementioned sensitivity of the LST to short-term synoptic variations, it seems that the approach is, therefore, not entirely robust to these issues. It is further important to mention that the generated T-S slope products, and by far most notably the LST trend image, exhibit a clear pattern of stripes that can be attributed to Landsat-7′s SLC-off error (see Figure 8 and Figure 9). For the LST image, slope coefficients may vary greatly even in adjacent homogenous areas, whereas the T-S trends of the NDVI and NDWI exhibit a far more stable picture. Due to this problem, the inherent variance of the LST and the comparably few observations, LST trends cannot solely rely on the trend results, and it seems highly advisable to investigate long and preferably dense time series to compensate for the apparently random information.
Nevertheless, as the correlation and regression results have revealed that an increase in mean greenness and especially moistness is associated with cooler LSTs, it can be inferred that the large-scale greening and wetting trends observed in the MDR will likely alter the surface thermal regime as indicated by the T-S slopes of the LST time series. The general Arctic warming trends, which are strongest during autumn/winter and lower in spring/summer [25,33], may, therefore, be expressed by cooler summer LSTs, as vegetation cover and moisture supply increases in the MDR.

5. Conclusions

This study investigated Landsat derived summer LST and multispectral indices between 1985 and 2018 and presents an overview of the mean summer LST, NDVI, and NDWI for the Arctic Mackenzie Delta Region, Northern Canada. The approach comprised the implementation of a Single-Channel algorithm within the GEE, which allows for a spatially flexible retrieval of LSTs over large areas and the calculation of long-term means to characterize the surface thermal regime at high spatial detail. The python script is available on GitHub (https://github.com/leonsnill/lst_landsat) and may serve as a valuable tool to derive statistical temporal metrics.
The correlation and regression analyses revealed that the predominant factors influencing the mean LST are the moisture conditions of the landscape that, in turn, are governed by the prevalent vegetation state and the topographic situation. The influence of vegetation and moisture properties on the LST is higher in terrain with subtle topography where differences in insolation are negligible, and the mean summer LST is best correlated with the mean TCW and the NDWI. Most of the study area’s land surface has been subject to a large-scale wetting and greening between 1985 and 2018. The positive trends of NDVI and NDWI are usually accompanied by negative trends of the summer LST. This inverse relationship is the result of intense shrubification and enhanced moisture supply that cools the surface via higher evapotranspiration rates with increased latent heat fluxes. Wildfire disturbance locally accelerates this process in the long-term as demonstrated for the 1968 wildfire area, which could be attributed to the intensive shrubification by deciduous species and the expansion of grasses that lead to significantly cooler LSTs than the undisturbed surroundings. The findings, therefore, depict the temporal dynamics and long-term results of the alteration of the surface, and it can be inferred that the large-scale greening and wetting trends observed in the MDR are cooling the mean summer LSTs. A multi-sensor approach could close the temporal gap between LST observations at high spatial detail, yet with a relatively low temporal resolution, and allow for absolute quantification of the alteration processes.

Author Contributions

Conceptualization, L.N. and T.U.; methodology, L.N. and T.U; software, L.N.; validation, L.N., T.U., C.K., J.S.-W. and R.B.; formal analysis, L.N. and T.U.; investigation, L.N. and T.U.; resources L.N., T.U., and R.B.; data curation, L.N. and T.U.; writing—original draft preparation, L.N. and T.U.; writing—review and editing, L.N., T.U., C.K., J.S.-W. and R.B.; visualization, L.N. and T.U.; supervision, T.U. and R.B.; project administration, T.U.; funding acquisition T.U., C.K., J.S.-W. and R.B.

Funding

This research was funded by the German Research Foundation (Deutsche Forschungsgemeinschaft—DFG) within the project “Multi-Scale Characterization of Polar Permafrost Landscapes by Airborne and Satellite Remote Sensing and In-Situ Geophysical Measurements” under grant number 329721376.

Acknowledgments

The Digital Elevation Model of the TanDEM-X Mission is shown under the permission of the German Aerospace Center (DLR), Germany—©DLR, 2016. The data were requested via the proposal “IDEM_HYDRO0182” and the TanDEM-X Science proposal “DEM_OTHER1323”. We thank Dr. Katharine Thomas for proofreading of the manuscript.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A

Table A1. Coefficients in matrix notation used to calculate the AFs in combination with water vapor contents. Coefficients are taken from [57,58].
Table A1. Coefficients in matrix notation used to calculate the AFs in combination with water vapor contents. Coefficients are taken from [57,58].
SensorCoefficients
Landsat-5 TM [ 0.07518 0.00492 1.03189 0.59600 1.22554 0.08104 0.02767 1.43740 0.25844 ]
Landsat-7 ETM+ [ 0.06518 0.00683 1.02717 0.53003 1.25866 0.10490 0.01965 1.36947 0.24310 ]
Landsat-8 OLI-TIRS [ 0.04019 0.02916 1.01523 0.38333 1.50294 0.20324 0.00918 1.36072 0.27514 ]

Appendix B

Table A2. Overview of the Landsat features derived in GEE and topography-related environmental factors derived from the TanDEM-X Digital Elevation Model.
Table A2. Overview of the Landsat features derived in GEE and topography-related environmental factors derived from the TanDEM-X Digital Elevation Model.
FeatureDescriptionYear/PeriodSource
LST_meanMean of Land Surface Temperature 1985−2018Landsat
LST_stdDevStandard Deviation of Land Surface Temperature 1985−2018Landsat
LST_tsTheil-Sen trend of Summer Land Surface Temperature 1985−2018Landsat
TCg_meanMean of Tasseled Cap Greenness 1985−2018Landsat
TCg_stdDevStandard Deviation of Tasseled Cap Greenness 1985−2018Landsat
TCg_tsTheil-Sen trend of Tasseled Cap Greenness 1985−2018Landsat
TCb_meanMean of Tasseled Cap Brightness 1985−2018Landsat
TCb_stdDevStandard Deviation of Tasseled Cap Brightness 1985−2018Landsat
TCb_tsTheil-Sen trend of Tasseled Cap Brightness 1985−2018Landsat
TCw_meanMean of Tasseled Cap Wetness 1985−2018Landsat
TCw_stdDevStandard Deviation of Tasseled Cap Wetness 1985−2018Landsat
TCw_tsTheil-Sen trend of Tasseled Cap Wetness 1985−2018Landsat
NDVI_meanMean of NDVI 1985−2018Landsat
NDVI_stdDevStandard Deviation of NDVI1985−2018Landsat
NDVI_tsTheil-Sen trend of NDVI1985−2018Landsat
NDWI_meanMean of NDWI1985−2018Landsat
NDWI_stdDevStandard Deviation of NDWI 1985−2018Landsat
NDWI_tsTheil-Sen trend of NDWI1985−2018Landsat
DEMTerrain Elevation2011−2012TanDEM-X
TWITopographic Wetness Index2011−2012TanDEM-X
FlowAccFlow Accumulation of Multi-Flow-Direction Approach 2011−2012TanDEM-X
InsolationPotential Annual Solar Insolation2011−2012TanDEM-X
STAGETerrain Exposition // Northness // Transformed Aspect 2011−2012TanDEM-X
WaterDistEuclidean Distance to WaterbodyVector Data

Appendix C

The tables in Appendix C summarize the multiple linear regression outputs as obtained from the statsmodels module in python. Statistics shown are the estimated coefficients (coef) and additionally obtained standardized coefficients (z-score coef), the standard errors (std err) of the coef’s, the t-statistic value (t) and associated p-value (P>|t|), the 95% confidence interval’s lower and upper values ([0.025 and 0.975], respectively), the Variance Inflation Factor (VIF) of each parameter, the explained variance (R2), the root-mean-square error (RMSE) and the Bayesian information criterion (BIC).
Table A3. Multiple linear regression output as obtained for the entire study area.
Table A3. Multiple linear regression output as obtained for the entire study area.
Study Area–Mackenzie Delta Region (Figure 6a)R2: 0.307RMSE: 1.730BIC: 2.351e+08
coefz-score coefstd errtP>|t|[0.0250.975]VIF
Intercept19.6059−2.234e−100.0013.23e+040.00019.60519.607
TCW_mean−23.4571−0.39140.007−3351.7130.000−23.471−23.4431.175
STAGE−41.6145−0.26870.004−2301.0840.000−9.071−9.0561.175
Table A4. Multiple linear regression output as obtained for the subregion “East Channel”.
Table A4. Multiple linear regression output as obtained for the subregion “East Channel”.
Subregion 1–East Channel (Figure 6b)R2: 0.726RMSE: 0.814BIC: 1.464e+06
coefz-score coefstd errtP>|t|[0.0250.975]VIF
Intercept17.5934−7.582e−080.0082251.6050.00017.57817.609
TCG_mean3.74840.06410.04583.6620.0003.6613.8361.294
TCW_mean−41.8219−0.6500.048−869.0130.000−41.916−41.7281.233
DEM0.00860.07840.00078.7540.0000.0080.0092.179
TWI−0.1850−0.22420.001−229.0990.000−0.187−0.1832.107
WaterDist0.00190.17877.77e−06243.9610.0000.0020.0021.180
Table A5. Multiple linear regression output as obtained for the subregion “Inuvik”.
Table A5. Multiple linear regression output as obtained for the subregion “Inuvik”.
Subregion 2–Inuvik (Figure 6c)R2: 0.629RMSE: 0.834BIC: 2.069e+06
coefz-score coefstd errtP>|t|[0.0250.975]VIF
Intercept21.5127−2.584e−150.0131693.9690.00021.48821.538
TCW_mean−52.6900−0.68550.080−660.3020.000−52.84−52.5342.426
NDVI_mean−4.5934−0.23870.017−277.9570.000−4.626−4.5611.659
TCB_mean−2.6552−0.05710.044−60.3000.000−2.741−2.5692.020
Table A6. Multiple linear regression output as obtained for the subregion “Delta”.
Table A6. Multiple linear regression output as obtained for the subregion “Delta”.
Subregion 3–Delta (Figure 6d)R2: 0.743RMSE: 0.511BIC: 7.426e+05
coefz-score coefstd errtP>|t|[0.0250.975]VIF
Intercept20.48626.883e−080.0073084.9680.00020.47320.499
NDWI_mean−10.3214−0.68460.015−690.6750.000−10.351−10.2921.894
TCW_mean−9.1922−0.16870.058−159.7250.000−9.305−9.0792.151
TWI0.08680.12080.001148.9720.0000.0860.0881.268
Table A7. Multiple linear regression output as obtained for the subregion “Richardson Mountains”.
Table A7. Multiple linear regression output as obtained for the subregion “Richardson Mountains”.
Subregion 4–Richardson Mountains (Figure 6e)R2: 0.479RMSE: 1.614BIC: 3.785e+06
coefz-score coefstd errtP>|t|[0.0250.975]VIF
Intercept18.5392−5.431e−090.0151220.1510.00018.50918.569
STAGE−11.9624−0.50220.024−506.3740.000−12.009−11.9161.882
TCB_mean10.57050.23940.039273.1890.00010.49510.6461.469
TCW_mean−2.6404−0.04470.051−51.6850.000−2.741−2.5401.429

References

  1. Climate Change 2013—The Physical Science Basis: Working Group I Contribution to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change; Intergovernmental Panel on Climate Change (Ed.) Cambridge University Press: Cambridge, UK, 2014; ISBN 978-1-107-41532-4. [Google Scholar]
  2. Nitze, I.; Grosse, G. Detection of landscape dynamics in the Arctic Lena Delta with temporally dense Landsat time-series stacks. Remote Sens. Environ. 2016, 181, 27–41. [Google Scholar] [CrossRef]
  3. Smith, S.L.; Romanovsky, V.E.; Lewkowicz, A.G.; Burn, C.R.; Allard, M.; Clow, G.D.; Yoshikawa, K.; Throop, J. Thermal state of permafrost in North America: A contribution to the international polar year. Permafr. Periglac. Process. 2010, 21, 117–135. [Google Scholar] [CrossRef]
  4. Lantz, T.C.; Kokelj, S.V. Increasing rates of retrogressive thaw slump activity in the Mackenzie Delta region, N.W.T., Canada. Geophys. Res. Lett. 2008, 35. [Google Scholar] [CrossRef]
  5. Fraser, R.; Olthof, I.; Carrière, M.; Deschamps, A.; Pouliot, D. A method for trend-based change analysis in Arctic tundra using the 25-year Landsat archive. Polar Rec. 2012, 48, 83–93. [Google Scholar] [CrossRef]
  6. Marsh, P.; Russell, M.; Pohl, S.; Haywood, H.; Onclin, C. Changes in thaw lake drainage in the Western Canadian Arctic from 1950 to 2000. Hydrol. Process. 2009, 23, 145–158. [Google Scholar] [CrossRef]
  7. Schuur, E. a. G.; McGuire, A.D.; Schädel, C.; Grosse, G.; Harden, J.W.; Hayes, D.J.; Hugelius, G.; Koven, C.D.; Kuhry, P.; Lawrence, D.M.; et al. Climate change and the permafrost carbon feedback. Nature 2015, 520, 171–179. [Google Scholar] [CrossRef]
  8. Hinzman, L.D.; Bettez, N.D.; Bolton, W.R.; Chapin, F.S.; Dyurgerov, M.B.; Fastie, C.L.; Griffith, B.; Hollister, R.D.; Hope, A.; Huntington, H.P.; et al. Evidence and Implications of Recent Climate Change in Northern Alaska and Other Arctic Regions. Clim. Chang. 2005, 72, 251–298. [Google Scholar] [CrossRef]
  9. Ju, J.; Masek, J.G. The vegetation greenness trend in Canada and US Alaska from 1984–2012 Landsat data. Remote Sens. Environ. 2016, 176, 1–16. [Google Scholar] [CrossRef]
  10. Fraser, R.H.; Lantz, T.C.; Olthof, I.; Kokelj, S.V.; Sims, R.A. Warming-Induced Shrub Expansion and Lichen Decline in the Western Canadian Arctic. Ecosystems 2014, 17, 1151–1168. [Google Scholar] [CrossRef]
  11. Chapin, F.S.; Sturm, M.; Serreze, M.C.; McFadden, J.P.; Key, J.R.; Lloyd, A.H.; McGuire, A.D.; Rupp, T.S.; Lynch, A.H.; Schimel, J.P.; et al. Role of Land-Surface Changes in Arctic Summer Warming. Science 2005, 310, 657–660. [Google Scholar] [CrossRef]
  12. Sturm, M.; Racine, C.; Tape, K. Increasing shrub abundance in the Arctic. Nature 2001, 411, 546–547. [Google Scholar] [CrossRef] [PubMed]
  13. Muster, S.; Langer, M.; Abnizova, A.; Young, K.L.; Boike, J. Spatio-temporal sensitivity of MODIS land surface temperature anomalies indicates high potential for large-scale land cover change detection in Arctic permafrost landscapes. Remote Sens. Environ. 2015, 168, 1–12. [Google Scholar] [CrossRef] [Green Version]
  14. Comiso, J.C. Warming Trends in the Arctic from Clear Sky Satellite Observations. J. Clim. 2003, 16, 3498–3510. [Google Scholar] [CrossRef] [Green Version]
  15. Comiso, J.C. Arctic warming signals from satellite observations. Weather 2006, 61, 70–76. [Google Scholar] [CrossRef] [Green Version]
  16. Langer, M.; Westermann, S.; Boike, J. Spatial and temporal variations of summer surface temperatures of wet polygonal tundra in Siberia - implications for MODIS LST based permafrost monitoring. Remote Sens. Environ. 2010, 114, 2059–2069. [Google Scholar] [CrossRef]
  17. Boike, J.; Grau, T.; Heim, B.; Günther, F.; Langer, M.; Muster, S.; Gouttevin, I.; Lange, S. Satellite-derived changes in the permafrost landscape of central Yakutia, 2000–2011: Wetting, drying, and fires. Glob. Planet. Chang. 2016, 139, 116–127. [Google Scholar] [CrossRef]
  18. Langer, M.; Westermann, S.; Heikenfeld, M.; Dorn, W.; Boike, J. Satellite-based modeling of permafrost temperatures in a tundra lowland landscape. Remote Sens. Environ. 2013, 135, 12–24. [Google Scholar] [CrossRef] [Green Version]
  19. Raynolds, M.K.; Comiso, J.C.; Walker, D.A.; Verbyla, D. Relationship between satellite-derived land surface temperatures, arctic vegetation types, and NDVI. Remote Sens. Environ. 2008, 112, 1884–1894. [Google Scholar] [CrossRef]
  20. Nguyen, T.-N.; Burn, C.R.; King, D.J.; Smith, S.L. Estimating the extent of near-surface permafrost using remote sensing, Mackenzie Delta, Northwest Territories. Permafr. Periglac. Process. 2009, 20, 141–153. [Google Scholar] [CrossRef]
  21. Lantz, T.C.; Kokelj, S.V.; Gergel, S.E.; Henry, G.H.R. Relative impacts of disturbance and temperature: Persistent changes in microenvironment and vegetation in retrogressive thaw slumps. Glob. Chang. Biol. 2009, 15, 1664–1675. [Google Scholar] [CrossRef]
  22. Lantz, T.C.; Gergel, S.E.; Henry, G.H.R. Response of green alder (Alnus viridis subsp. fruticosa) patch dynamics and plant community composition to fire and regional temperature in north-western Canada. J. Biogeogr. 2010, 37, 1597–1610. [Google Scholar] [CrossRef]
  23. Smith, S.L.; Burgess, M.M.; Riseborough, D.; Nixon, F.M. Recent trends from Canadian permafrost thermal monitoring network sites. Permafr. Periglac. Process. 2005, 16, 19–30. [Google Scholar] [CrossRef]
  24. Lantz, T.C.; Marsh, P.; Kokelj, S.V. Recent Shrub Proliferation in the Mackenzie Delta Uplands and Microclimatic Implications. Ecosystems 2013, 16, 47–59. [Google Scholar] [CrossRef]
  25. Burn, C.R.; Kokelj, S.V. The environment and permafrost of the Mackenzie Delta area. Permafr. Periglac. Process. 2009, 20, 83–105. [Google Scholar] [CrossRef]
  26. Fraser, R.H.; Olthof, I.; Kokelj, S.V.; Lantz, T.C.; Lacelle, D.; Brooker, A.; Wolfe, S.; Schwarz, S. Detecting Landscape Changes in High Latitude Environments Using Landsat Trend Analysis: 1. Visualization. Remote Sens. 2014, 6, 11533–11557. [Google Scholar] [CrossRef] [Green Version]
  27. Eugster, W.; Rouse, W.R.; Pielke Sr, R.A.; Mcfadden, J.P.; Baldocchi, D.D.; Kittel, T.G.F.; Chapin, F.S.; Liston, G.E.; Vidale, P.L.; Vaganov, E.; et al. Land-atmosphere energy exchange in Arctic tundra and boreal forest: Available data and feedbacks to climate. Glob. Chang. Biol. 2000, 6, 84–115. [Google Scholar] [CrossRef]
  28. Goulding, H.L.; Prowse, T.D.; Bonsal, B. Hydroclimatic controls on the occurrence of break-up and ice-jam flooding in the Mackenzie Delta, NWT, Canada. J. Hydrol. 2009, 379, 251–267. [Google Scholar] [CrossRef]
  29. MacKay, J.R. The Mackenzie Delta Area, N.W.T.; Queen’s Printer: Ottawa, ON, Canada, 1963. [Google Scholar]
  30. Lantz, T.C.; Gergel, S.E.; Kokelj, S.V. Spatial Heterogeneity in the Shrub Tundra Ecotone in the Mackenzie Delta Region, Northwest Territories: Implications for Arctic Environmental Change. Ecosystems 2010, 13, 194–204. [Google Scholar] [CrossRef]
  31. Kokelj, S.V.; Lantz, T.C.; Kanigan, J.; Smith, S.L.; Coutts, R. Origin and polycyclic behaviour of tundra thaw slumps, Mackenzie Delta region, Northwest Territories, Canada. Permafr. Periglac. Process. 2009, 20, 173–184. [Google Scholar] [CrossRef]
  32. Olthof, I.; Fraser, R.H. Detecting Landscape Changes in High Latitude Environments Using Landsat Trend Analysis: 2. Classification. Remote Sens. 2014, 6, 11558–11578. [Google Scholar] [CrossRef] [Green Version]
  33. Serreze, M.C.; Barry, R.G. Processes and impacts of Arctic amplification: A research synthesis. Glob. Planet. Chang. 2011, 77, 85–96. [Google Scholar] [CrossRef]
  34. Gillett, N.P.; Weaver, A.J.; Zwiers, F.W.; Flannigan, M.D. Detecting the effect of climate change on Canadian forest fires. Geophys. Res. Lett. 2004, 31. [Google Scholar] [CrossRef] [Green Version]
  35. Hu, F.S.; Higuera, P.E.; Walsh, J.E.; Chapman, W.L.; Duffy, P.A.; Brubaker, L.B.; Chipman, M.L. Tundra burning in Alaska: Linkages to climatic change and sea ice retreat. J. Geophys. Res. Biogeosci. 2010, 115. [Google Scholar] [CrossRef] [Green Version]
  36. McCoy, V.M.; Burn, C.R. Potential Alteration by Climate Change of the Forest-Fire Regime in the Boreal Forest of Central Yukon Territory; Arctic Institute of North America: Calgary, AB, Canada, 2010. [Google Scholar]
  37. Jones, B.M.; Grosse, G.; Arp, C.D.; Miller, E.; Liu, L.; Hayes, D.J.; Larsen, C.F. Recent Arctic tundra fire initiates widespread thermokarst development. Sci. Rep. 2015, 5. [Google Scholar] [CrossRef] [PubMed]
  38. Gorelick, N.; Hancher, M.; Dixon, M.; Ilyushchenko, S.; Thau, D.; Moore, R. Google Earth Engine: Planetary-scale geospatial analysis for everyone. Remote Sens. Environ. 2017, 202, 18–27. [Google Scholar] [CrossRef]
  39. Krieger, G.; Moreira, A.; Fiedler, H.; Hajnsek, I.; Werner, M.; Younis, M.; Zink, M. TanDEM-X: A Satellite Formation for High-Resolution SAR Interferometry. Ieee Trans. Geosci. Remote Sens. 2007, 45, 3317–3341. [Google Scholar] [CrossRef] [Green Version]
  40. Gruber, A.; Wessel, B.; Huber, M.; Roth, A. Operational TanDEM-X DEM calibration and first validation results. Isprs J. Photogramm. Remote Sens. 2012, 73, 39–49. [Google Scholar] [CrossRef]
  41. Deutsches Zentrum für Luft- und Raumfahrt (DLR) TanDEM-X Ground Segment DEM Products Specification Document 2016; EOC–Earth Observation Center: Weßling, Germany, 2016.
  42. Deutsches Zentrum für Luft- und Raumfahrt (DLR) TanDEM-X Science Service System 2018; EOC–Earth Observation Center: Weßling, Germany, 2018.
  43. Wilson, J.P.; Gallant, J.C. Terrain Analysis: Principles and Applications; Wiley: New York, NY, USA, 2000; ISBN 978-0-471-32188-0. [Google Scholar]
  44. Hofierka, J.; Súri, M.; Marecka, M. The solar radiation model for Open source GIS: Implementation and applications. In Proceedings of the Open source GIS-GRASS users conference 2002, Trento, Italy, 11–13 September 2002. [Google Scholar]
  45. Gruber, S.; Peckham, S. Chapter 7 Land-Surface Parameters and Objects in Hydrology. In Developments in Soil Science; Hengl, T., Reuter, H.I., Eds.; Geomorphometry; Elsevier: Amsterdam, The Netherlands, 2009; Volume 33, pp. 171–194. [Google Scholar]
  46. Rouse, J.W. Monitoring Vegetation Systems in the Great Plains with ERTS; NASA Technical Reports Server, 1974. [Google Scholar]
  47. Gao, B. NDWI—A normalized difference water index for remote sensing of vegetation liquid water from space. Remote Sens. Environ. 1996, 58, 257–266. [Google Scholar] [CrossRef]
  48. McFEETERS, S.K. The use of the Normalized Difference Water Index (NDWI) in the delineation of open water features. Int. J. Remote Sens. 1996, 17, 1425–1432. [Google Scholar] [CrossRef]
  49. Ji, L.; Zhang, L.; Wylie, B.K.; Rover, J. On the terminology of the spectral vegetation index (NIR − SWIR)/(NIR + SWIR). Int. J. Remote Sens. 2011, 32, 6901–6909. [Google Scholar] [CrossRef]
  50. Nitze, I.; Grosse, G.; Jones, B.M.; Romanovsky, V.E.; Boike, J. Remote sensing quantifies widespread abundance of permafrost region disturbances across the Arctic and Subarctic. Nat Commun 2018, 9. [Google Scholar] [CrossRef] [PubMed]
  51. Chander, G.; Markham, B. Revised Landsat-5 TM radiometric calibration procedures and postcalibration dynamic ranges. Ieee Trans. Geosci. Remote Sens. 2003, 41, 2674–2677. [Google Scholar] [CrossRef]
  52. Chander, G.; Markham, B.L.; Helder, D.L. Summary of current radiometric calibration coefficients for Landsat MSS, TM, ETM+, and EO-1 ALI sensors. Remote Sens. Environ. 2009, 113, 893–903. [Google Scholar] [CrossRef]
  53. Qin, Z.; Karnieli, A.; Berliner, P. A mono-window algorithm for retrieving land surface temperature from Landsat TM data and its application to the Israel-Egypt border region. Int. J. Remote Sens. 2001, 22, 3719–3746. [Google Scholar] [CrossRef]
  54. Jiménez-Muñoz, J.C.; Sobrino, J.A. A generalized single-channel method for retrieving land surface temperature from remote sensing data. J. Geophys. Res. Atmos. 2003, 108. [Google Scholar] [CrossRef]
  55. Montanaro, M.; Gerace, A.; Lunsford, A.; Reuter, D. Stray Light Artifacts in Imagery from the Landsat 8 Thermal Infrared Sensor. Remote Sens. 2014, 6, 10435–10456. [Google Scholar] [CrossRef] [Green Version]
  56. Loveland, T.R.; Irons, J.R. Landsat 8: The plans, the reality, and the legacy. Remote Sens. Environ. 2016, 185, 1–6. [Google Scholar] [CrossRef] [Green Version]
  57. Jimenez-Munoz, J.C.; Cristobal, J.; Sobrino, J.A.; Soria, G.; Ninyerola, M.; Pons, X.; Pons, X. Revision of the Single-Channel Algorithm for Land Surface Temperature Retrieval From Landsat Thermal-Infrared Data. Ieee Trans. Geosci. Remote Sens. 2009, 47, 339–349. [Google Scholar] [CrossRef]
  58. Jiménez-Muñoz, J.C.; Sobrino, J.A.; Skoković, D.; Mattar, C.; Cristóbal, J. Land Surface Temperature Retrieval Methods From Landsat-8 Thermal Infrared Sensor Data. Ieee Geosci. Remote Sens. Lett. 2014, 11, 1840–1843. [Google Scholar] [CrossRef]
  59. Kalnay, E.; Kanamitsu, M.; Kistler, R.; Collins, W.; Deaven, D.; Gandin, L.; Iredell, M.; Saha, S.; White, G.; Woollen, J.; et al. The NCEP/NCAR 40-Year Reanalysis Project. Bull. Amer. Meteor. Soc. 1996, 77, 437–472. [Google Scholar] [CrossRef]
  60. Li, Z.-L.; Tang, B.-H.; Wu, H.; Ren, H.; Yan, G.; Wan, Z.; Trigo, I.F.; Sobrino, J.A. Satellite-derived land surface temperature: Current status and perspectives. Remote Sens. Environ. 2013, 131, 14–37. [Google Scholar] [CrossRef] [Green Version]
  61. Rosas, J.; Houborg, R.; McCabe, M.F. Sensitivity of Landsat 8 Surface Temperature Estimates to Atmospheric Profile Data: A Study Using MODTRAN in Dryland Irrigated Systems. Remote Sens. 2017, 9, 988. [Google Scholar] [CrossRef]
  62. Sobrino, J.A.; Jimenez-Munoz, J.C.; Soria, G.; Romaguera, M.; Guanter, L.; Moreno, J.; Plaza, A.; Martinez, P. Land Surface Emissivity Retrieval From Different VNIR and TIR Sensors. Ieee Trans. Geosci. Remote Sens. 2008, 46, 316–327. [Google Scholar] [CrossRef]
  63. Li, Z.-L.; Wu, H.; Wang, N.; Qiu, S.; Sobrino, J.A.; Wan, Z.; Tang, B.-H.; Yan, G. Land surface emissivity retrieval from satellite data. Int. J. Remote Sens. 2013, 34, 3084–3127. [Google Scholar] [CrossRef]
  64. Carlson, T.N.; Ripley, D.A. On the relation between NDVI, fractional vegetation cover, and leaf area index. Remote Sens. Environ. 1997, 62, 241–252. [Google Scholar] [CrossRef]
  65. Dash, P.; Göttsche, F.-M.; Olesen, F.-S.; Fischer, H. Separating surface emissivity and temperature using two-channel spectral indices and emissivity composites and comparison with a vegetation fraction method. Remote Sens. Environ. 2005, 96, 1–17. [Google Scholar] [CrossRef]
  66. Sobrino, J.A.; Jiménez-Muñoz, J.C.; Paolini, L. Land surface temperature retrieval from LANDSAT TM 5. Remote Sens. Environ. 2004, 90, 434–440. [Google Scholar] [CrossRef]
  67. Wang, F.; Qin, Z.; Song, C.; Tu, L.; Karnieli, A.; Zhao, S. An Improved Mono-Window Algorithm for Land Surface Temperature Retrieval from Landsat 8 Thermal Infrared Sensor Data. Remote Sens. 2015, 7, 4268–4289. [Google Scholar] [CrossRef] [Green Version]
  68. Brooker, A.; Fraser, R.H.; Olthof, I.; Kokelj, S.V.; Lacelle, D. Mapping the Activity and Evolution of Retrogressive Thaw Slumps by Tasselled Cap Trend Analysis of a Landsat Satellite Image Stack. Permafr. Periglac. Process. 2014, 25, 243–256. [Google Scholar] [CrossRef]
  69. Nitze, I.; Grosse, G.; Jones, B.M.; Arp, C.D.; Ulrich, M.; Fedorov, A.; Veremeeva, A. Landsat-Based Trend Analysis of Lake Dynamics across Northern Permafrost Regions. Remote Sens. 2017, 9, 640. [Google Scholar] [CrossRef]
  70. Sen, P.K. Estimates of the Regression Coefficient Based on Kendall’s Tau. J. Am. Stat. Assoc. 1968, 63, 1379–1389. [Google Scholar] [CrossRef]
  71. Bhatt, U.S.; Walker, D.A.; Raynolds, M.K.; Comiso, J.C.; Epstein, H.E.; Jia, G.; Gens, R.; Pinzon, J.E.; Tucker, C.J.; Tweedie, C.E.; et al. Circumpolar Arctic Tundra Vegetation Change Is Linked to Sea Ice Decline. Earth Interact. 2010, 14, 1–20. [Google Scholar] [CrossRef] [Green Version]
  72. Mackay, J.R. Active Layer Changes (1968 to 1993) following the Forest-Tundra Fire near Inuvik, N.W.T., Canada. Arct. Alp. Res. 1995, 27, 323–336. [Google Scholar] [CrossRef]
  73. Jackson, T.J.; Chen, D.; Cosh, M.; Li, F.; Anderson, M.; Walthall, C.; Doriaswamy, P.; Hunt, E.R. Vegetation water content mapping using Landsat data derived normalized difference water index for corn and soybeans. Remote Sens. Environ. 2004, 92, 475–482. [Google Scholar] [CrossRef]
  74. Bhatt, U.S.; Walker, D.A.; Raynolds, M.K.; Bieniek, P.A.; Epstein, H.E.; Comiso, J.C.; Pinzon, J.E.; Tucker, C.J.; Polyakov, I.V. Recent Declines in Warming and Vegetation Greening Trends over Pan-Arctic Tundra. Remote Sens. 2013, 5, 4229–4254. [Google Scholar] [CrossRef] [Green Version]
  75. Markus, T.; Stroeve, J.C.; Miller, J. Recent changes in Arctic sea ice melt onset, freezeup, and melt season length. J. Geophys. Res. Ocean. 2009, 114. [Google Scholar] [CrossRef]
  76. Wein, R.W. Forest Fires and Northern Communities: Lessons from the 1968 Inuvik fir 2002; University of Alberta: Inuvik, NWT, Canada, 2002. [Google Scholar]
Figure 1. Mackenzie Delta Region: (a) Landsat-8 RGB true-color image composite of the study area and (b) overview. Subregions are highlighted by red polygons in (a), and photographs) of the subregions are numbered from 1–4 (bottom). Subfigure (c) displays the TanDEM-X digital elevation model (© DLR 2016) superimposed by the water-mask (black). Photography by T. Ullmann 2012, 2013, 2018.
Figure 1. Mackenzie Delta Region: (a) Landsat-8 RGB true-color image composite of the study area and (b) overview. Subregions are highlighted by red polygons in (a), and photographs) of the subregions are numbered from 1–4 (bottom). Subfigure (c) displays the TanDEM-X digital elevation model (© DLR 2016) superimposed by the water-mask (black). Photography by T. Ullmann 2012, 2013, 2018.
Remotesensing 11 02329 g001
Figure 2. Landsat observations: (a) pixel-based number of cloud-free summer (July/August) observations between 1985 and 2018; (b) number of scenes per month. Subplots from Figure 1 are highlighted by the black polygons.
Figure 2. Landsat observations: (a) pixel-based number of cloud-free summer (July/August) observations between 1985 and 2018; (b) number of scenes per month. Subplots from Figure 1 are highlighted by the black polygons.
Remotesensing 11 02329 g002
Figure 3. Flowchart of the processing steps, including pre-processing the data, retrieving Land Surface Temperature (LST), and calculating the statistical parameters in order to conduct the correlation, regression, and trend analysis.
Figure 3. Flowchart of the processing steps, including pre-processing the data, retrieving Land Surface Temperature (LST), and calculating the statistical parameters in order to conduct the correlation, regression, and trend analysis.
Remotesensing 11 02329 g003
Figure 4. Pixel based means (ac) and standard deviations (df) of summer (i.e., July, August) Land Surface Temperature (LST), Normalized Difference Vegetation Index (NDVI), and Normalized Difference Water Index (NDWI) for the period from 1985 to 2018. The red rectangular boxes highlight the subregions, and the white polygons indicate the extent of former wildfires.
Figure 4. Pixel based means (ac) and standard deviations (df) of summer (i.e., July, August) Land Surface Temperature (LST), Normalized Difference Vegetation Index (NDVI), and Normalized Difference Water Index (NDWI) for the period from 1985 to 2018. The red rectangular boxes highlight the subregions, and the white polygons indicate the extent of former wildfires.
Remotesensing 11 02329 g004
Figure 5. Subsets of the correlation analyses as linear Pearson Correlation Coefficients (R) for the entire study area, and the subplots “East Channel”, “Inuvik”, “Delta”, “Richardson Mountains”. A tabular summary of all investigated parameters and their abbreviations is provided in Appendix B (Table A2).
Figure 5. Subsets of the correlation analyses as linear Pearson Correlation Coefficients (R) for the entire study area, and the subplots “East Channel”, “Inuvik”, “Delta”, “Richardson Mountains”. A tabular summary of all investigated parameters and their abbreviations is provided in Appendix B (Table A2).
Remotesensing 11 02329 g005
Figure 6. Predicted-observed density plots of the multiple linear regressions for the entire study area (a) and the four subregions (be). Detailed information on each model is provided in Appendix C (Table A3, Table A4, Table A5, Table A6 and Table A7).
Figure 6. Predicted-observed density plots of the multiple linear regressions for the entire study area (a) and the four subregions (be). Detailed information on each model is provided in Appendix C (Table A3, Table A4, Table A5, Table A6 and Table A7).
Remotesensing 11 02329 g006
Figure 7. Theil-Sen trend slopes of (a) summer LST, (b) NDVI, and (c) NDWI. The boundaries of areas affected by former wildfires are highlighted in white, whereas the red boxes highlight the subregions shown in Figure 8 and Figure 9.
Figure 7. Theil-Sen trend slopes of (a) summer LST, (b) NDVI, and (c) NDWI. The boundaries of areas affected by former wildfires are highlighted in white, whereas the red boxes highlight the subregions shown in Figure 8 and Figure 9.
Remotesensing 11 02329 g007
Figure 8. Local temporal dynamics for the “East Channel” subregion: true color RGB Landsat image (a), Theil-Sen slopes of the LST (b), NDVI (c), and NDWI (d), as well as Theil-Sen regression lines of selected locations (eg). Each plot location represents the spatial average of 90 × 90 m. The trend lines are complemented by the lower and upper limit of the 95% confidence interval of the slope coefficient.
Figure 8. Local temporal dynamics for the “East Channel” subregion: true color RGB Landsat image (a), Theil-Sen slopes of the LST (b), NDVI (c), and NDWI (d), as well as Theil-Sen regression lines of selected locations (eg). Each plot location represents the spatial average of 90 × 90 m. The trend lines are complemented by the lower and upper limit of the 95% confidence interval of the slope coefficient.
Remotesensing 11 02329 g008
Figure 9. Local temporal dynamics for the “Inuvik” subregion: true color RGB Landsat image (a) Theil-Sen slopes of the LST (b) NDVI (c) and NDWI (d) as well as Theil-Sen regression lines of selected locations (eg). Each plot location represents the spatial average of 90 × 90 m. The trend lines are complemented by the lower and upper limit of the 95% confidence interval of the slope coefficient. The dashed line in (g) indicates the 2012 wildfire disturbance.
Figure 9. Local temporal dynamics for the “Inuvik” subregion: true color RGB Landsat image (a) Theil-Sen slopes of the LST (b) NDVI (c) and NDWI (d) as well as Theil-Sen regression lines of selected locations (eg). Each plot location represents the spatial average of 90 × 90 m. The trend lines are complemented by the lower and upper limit of the 95% confidence interval of the slope coefficient. The dashed line in (g) indicates the 2012 wildfire disturbance.
Remotesensing 11 02329 g009

Share and Cite

MDPI and ACS Style

Nill, L.; Ullmann, T.; Kneisel, C.; Sobiech-Wolf, J.; Baumhauer, R. Assessing Spatiotemporal Variations of Landsat Land Surface Temperature and Multispectral Indices in the Arctic Mackenzie Delta Region between 1985 and 2018. Remote Sens. 2019, 11, 2329. https://doi.org/10.3390/rs11192329

AMA Style

Nill L, Ullmann T, Kneisel C, Sobiech-Wolf J, Baumhauer R. Assessing Spatiotemporal Variations of Landsat Land Surface Temperature and Multispectral Indices in the Arctic Mackenzie Delta Region between 1985 and 2018. Remote Sensing. 2019; 11(19):2329. https://doi.org/10.3390/rs11192329

Chicago/Turabian Style

Nill, Leon, Tobias Ullmann, Christof Kneisel, Jennifer Sobiech-Wolf, and Roland Baumhauer. 2019. "Assessing Spatiotemporal Variations of Landsat Land Surface Temperature and Multispectral Indices in the Arctic Mackenzie Delta Region between 1985 and 2018" Remote Sensing 11, no. 19: 2329. https://doi.org/10.3390/rs11192329

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop