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

: 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 surﬁcial properties and characterizes the land–atmosphere exchange of energy. Accordingly, this study analyses the inﬂuence of di ﬀ erent 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 Di ﬀ erence Vegetation Index (NDVI), Normalized Di ﬀ erence 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 inﬂuence 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 inﬂuenced by the topographic exposition as well as the prevalent moisture regime where higher evapotranspiration rates increase the latent heat ﬂux 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 di ﬀ erent regions of the MDR (R 2 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 inﬂuential factors and suggests accounting for a broad spectrum of di ﬀ erent 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 wildﬁres intensiﬁes the surface alterations locally and lead to signiﬁcantly cooler LSTs in the long-term compared to the undisturbed surroundings. and TCB, the only shows on mean variance,


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 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.

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 km 2 [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 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 Remote Sens. 2019, 11, 2329 5 of 26 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].

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). 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].

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 Remote Sens. 2019, 11, 2329 6 of 26 (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 2 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. 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). 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 parameterthe 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. 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). 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.

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]: 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]: 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)): where the sensor and band-specific coefficients C x that were used in this study are summarized in [2], and ρ x are the reflectance values in the corresponding parts of the electromagnetic spectrum.

Retrieval of Land Surface Temperature
For Landsat thermal infrared data with a certain channel width, one can obtain an effective at-satellite temperature BT sen (Equation (4)) based on an approximation of Planck's law as follows [51]: where K 1 and K 2 are band-specific thermal conversion constants provided with the metadata and L sen refers to the spectral radiance in W/(m 2 ·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)): 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]: where L sen is the spectral radiance at the sensor in W/m 2 ·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: where BT sen 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): The coefficients c ij 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 (SNDVI THM ) 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 NDVI v ), as well as bare soil (ε s and NDVI s ). The pixel emissivity values may then be derived accordingly [62] (Equation (10)): Remote Sens. 2019, 11, 2329 9 of 26 with P V being the fraction of vegetation cover [64] (Equation (11)): NDVI s and NDVI v 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)).

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]. 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 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 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".

Mean Summer Land Surface Temperature, NDVI, and NDWI
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.

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 (Tables A3-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.
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 R 2 of 0.73 and an RMSE of 0.81 °C. By comparing the standardized coefficient estimates (Table C1-C5), 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. 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).
Remote Sens. 2019, 10, x FOR PEER REVIEW 13 of 25 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 B1).

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 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 R 2 of 0.73 and an RMSE of 0.81 • C. By comparing the standardized coefficient estimates (Tables A3-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.

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. 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.

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.

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

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.

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].

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.

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 Figures 8 and 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.

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 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. Table A1. Coefficients in matrix notation used to calculate the AFs in combination with water vapor contents. Coefficients are taken from [57,58].

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 (R 2 ), the root-mean-square error (RMSE) and the Bayesian information criterion (BIC).