Towards physical habitat characterisation in the Antarctic S ø r Rondane Mountains using satellite remote sensing

Ice-free areas occupy less than 0.25% of the Antarctic surface, and mainly occur along coastlines, or as inland nunataks protruding from the extensive ice sheet. Their extreme environment and geographical isolation have contributed to the evolution of highly adapted, and largely endemic terrestrial biota. Physical habitat mapping is important to identify the main drivers of spatial variation in soil biodiversity, to predict its response to climate change and to help conservation planning. In this paper we retrieved remotely sensed Land Surface Temperature (LST) and Digital Surface Models (DSM) for the S ø r Rondane Mountains in East Antarctica from respectively the Thermal InfraRed Sensor (TIRS) on Landsat 8 and the Pl ´ eiades constellation of high resolution optical imagers. Satellite data were combined with ground truth temperature and elevation measurements with the aim to assess the performance of these remotely sensed data. Over a 2 year period, satellite derived LST corresponded to in situ temperature with Mean Average

Ice-free areas occupy less than 0.25% of the Antarctic surface, and mainly occur along coastlines, or as inland nunataks protruding from the extensive ice sheet.Their extreme environment and geographical isolation have contributed to the evolution of highly adapted, and largely endemic terrestrial biota.Physical habitat mapping is important to identify the main drivers of spatial variation in soil biodiversity, to predict its response to climate change and to help conservation planning.In this paper we retrieved remotely sensed Land Surface Temperature (LST) and Digital Surface Models (DSM) for the Sør Rondane Mountains in East Antarctica from respectively the Thermal InfraRed Sensor (TIRS) on Landsat 8 and the Pléiades constellation of high resolution optical imagers.Satellite data were combined with ground truth temperature and elevation measurements with the aim to assess the performance of these remotely sensed data.Over a 2 year period, satellite derived LST corresponded to in situ temperature with Mean Average Difference (MAD) of − 2.5 K, and Root Mean Squared Difference (RMSD) of 6.3 K. Lower biases were observed for periods when the data loggers were frozen or snow covered (MAD − 1.8K) compared to intervals where the devices were not frozen (MAD − 4 K), with larger scatter (RMSD 7 K) being observed during the latter.These larger MAD and RMSD are caused by the differential heating of the top of the gravel/rocks as observed by the satellite compared to their underside, where the loggers were installed, and are also impacted by the time step of the in situ logging (3 h) and the non-linear heating of the surface.In addition, for the different study sites different biases were observed as a result of the spatial resolution of the TIRS, depending on the composition, structure, geomorphology, and the surroundings of the logger position.Sites on a narrow rock outcrop (e.g. the Perlebandet and Utsteinen nunataks) show a negative bias due to the surrounding ice fields decreasing the satellite pixel average temperature compared to the in situ measured temperatures.For sites with a more extensive rock and gravel composition further away from the ice sheet (Dry and Yuboku Valleys), a positive bias was found as a result of the temperature differential between the exposed top and covered bottom of the rocks and gravel.The DSM derived from Pléiades (without bias correction) showed in general a bias of − 6 m (MAD) and scatter of 9 m (RMSD) when compared to the Reference Elevation Model of Antarctica and IceSAT-1 LIDAR data.Lowest errors were also found here for the extensive Dry and Yuboku Valley sites (respectively − 5 and 3 m MAD, 7 and 4 m RMSD).Correspondence with in situ logged GPS positions was slightly worse, with a positive bias of 12 m (MAD) and scatter of 15 m (RMSD) across the sites.An analysis on the basis of these LST and DSM datasets, show a clear separability of sites by their average temperature and solar exposition (in terms of slope and aspect) at the time of Landsat overpass.Due to the Landsat overpass times, the LST datasets are however strongly biased, and further work is needed for a better understanding of the light and temperature climate in the ice-free regions.This could for example be achieved by coupling the DSM to a solar irradiance and terrain shadowing model.

Introduction
Less than a quarter percent of the Antarctic continent consists of icefree areas (Burton-Johnson et al., 2016), which are confined to coastal regions and nunataks that protrude from the ice sheet (Cary et al., 2010).Environmental conditions in these ice-free areas are among the most extreme on Earth, and life is therefore dominated by microorganisms (Pushkareva et al., 2018) that show a high degree of endemism (Vyverman et al., 2010;Fernandez-Carazo et al., 2012;Verleyen et al., 2021).Microbial community structure is strongly dependent on received solar radiation (Convey et al., 2014), and available soil moisture (Obbels et al., 2016;Tytgat et al., 2016;Ray et al., 2020;Ortiz et al., 2020) during the brief austral summer.Soil moisture is in turn is influenced by soil temperature, upon which the metabolic activity also depends.The communities are structured by these microclimatic conditions as a result of the (micro)topography of habitats.The distribution of organisms may vary significantly at short spatial scales, and can include a vertical component, e.g.organisms are located not on, but just below the surface (Davey and Clarke, 1991), likely to avoid drier surface conditions (Davey et al., 1992).The slope of the terrain influences on the one hand the received solar irradiance which is maximal for a surface orthogonal to the incoming light, and on the other hand affects the long-term stability of the habitat.Similarly, sun exposure is a function of the orientation of the slope relative to the path of the sun through the sky.Combined with the regional topography and wind patterns that influence the presence of snow, these factors regulate soil moisture and temperature.Drier areas are characterised by sparse microbial soil crusts and biofilms that grow slowly and are very sensitive to terrain and climate conditions (Tytgat et al., 2016).Comparatively, wet areas are oases where macroscopic cryptogams such as mosses and lichens as well as microbial mats consisting of filamentous cyanobacteria and chlorophytes may develop (Obbels et al., 2016).These areas with visible life are found along ephemeral streams, under melting snowbanks, around the edges of lakes or ponds, and inside cryoconite holes (Priscu et al., 1998;Paerl and Priscu, 1998;McKnight et al., 1999;Niederberger et al., 2012;Lutz et al., 2019;Sohm et al., 2020;Weisleitner et al., 2020).
The surface temperature in polar regions varies greatly at short temporal and spatial scales, and needs to be robustly characterised in order to understand its effect on biota, before even attempting analysing impacts of e.g.climate change (Davey et al., 1992;Convey et al., 2018).The in situ characterisation of the exposure, temperature and humidity of these ice-free areas is difficult, expensive, and time consuming.The use of satellite data may provide a significant advantage for large scale characterisation of these areas.Mesoscale model outputs have been used for broad scale classification of the continent's ice-free areas, e.g. in 200 × 200 km cells by Terauds et al. (2012), but due to the presence of very strong spatial heterogeneity (Convey et al., 2014), the use of high resolution satellite remote sensing, e.g. from thermal imagers on Landsat or Aqua/Terra platforms, may be more appropriate for the characterisation of local temperature gradients.As temperature drives the liquid water cycle and metabolic rates, with both direct and indirect effects on microbial community structure and growth, satellite derived Land Surface Temperatures (LST) may be of use for differentiating between areas with different types of soil microbiomes.
Many authors have used the Landsat series of satellites to derive Land (LST) and Water Surface Temperatures (WST), e.g. for terrestrial (Yuan and Bauer, 2007), and coastal monitoring (Thomas et al., 2002), but also for studying the Antarctic ice shelves and sea ice (Scambos et al., 2000(Scambos et al., , 2004)).Orheim and Lucchitta (1988) used Landsat 5 imagery to estimate surface temperatures of the ice sheet in Dronning Maud Land.More recently, the lowest surface temperatures on earth were found on the snow surfaces between the east Antarctic Dome A to F using Landsat 8/TIRS imagery (Scambos et al., 2018).Thermal data from Landsat 8/TIRS and MODIS were used for modeling melt rates of Antarctic glaciers (Brabyn and Stichbury, 2020).Optical data from Landsat have been used for the construction of the Landsat Image Mosaic of Antarctica (LIMA, Bindschadler et al., 2008), and for studies of ice flow rates (Hambrey and Dowdeswell, 1994;Scambos et al., 2004;Fahnestock et al., 2016) and snow characteristics (Bourdelles and Fily, 1993).Rock outcrops have been mapped at a continental scale using Landsat 8 imagery (Burton-Johnson et al., 2016), yet applications of Landsat thermal data in the ice-free areas of the Antarctic continent are scarce.One study found that the LST derived from Landsat 7/ETM + compared to in situ logged temperatures with a 3.8K Root Mean Squared Error for study sites in the McMurdo Dry Valleys (Brabyn et al., 2014).
There are a wide range of applications of Digital Elevation Models or Digital Surface Models (DSM), including for example, the development of watershed models (Band, 1986), monitoring of erosion (Martınez-Casasnovas et al., 2002) and earthquake impacts (Van Puymbroeck et al., 2000;Turker and Cetinkaya, 2005), flood risk mapping and coastline modelling (Bates and De Roo, 2000;Demirkesen et al., 2007;Kulp and Strauss, 2018), and studies of forest canopy dynamics (Miller et al., 2000;Henbo et al., 2006).Applications in polar regions are of particular interest, for example to document ice sheet extent and mass balance (Rignot and Thomas, 2002;Shepherd et al., 2012), but also for modeling of available light in mountainous areas (Hansen et al., 2002) or in the McMurdo Dry Valleys (Dana et al., 1998;Acosta et al., 2020).The Bedmap1 (Lythe and Vaughan, 2001) and Bedmap2 (Fretwell et al., 2013) datasets cover the Antarctic continent with a 1 km grid, and include a digital elevation model, but their focus lies on subglacial bed elevation and ice-sheet thickness.The authors stress that the ice-sheet thickness was originally gridded at 5 km, and even at that resolution, many grid cells (66%) have no reference data (Lythe and Vaughan, 2001;Fretwell et al., 2013), illustrating the challenge of continent-wide mapping efforts.The Bedmap2 elevation data are based on the Liu et al. (1999) Antarctica elevation model, with an assigned uncertainty of 30 m in general, and 130 m in mountainous areas (Liu et al., 1999;Fretwell et al., 2013).The Liu et al. (1999) model has variable horizontal resolution of 200 m in the mountains, 400 m in the coastal areas, and 5 km in the interior, and combines a variety of elevation sources, including cartographic data, satellite altimetry, and GPS traverse data.Mapping of terrain elevation from optical satellite imagery generally uses the parallax information contained in stereo acquisitions made under different viewing conditions, ideally under the same illumination conditions.Several commercially exploited satellites are able to provide high resolution imagery for this application, e.g. the WorldView and Pléiades series of satellites.Challenges for high latitude and mountainous areas with steep relief are the lack of contrast or the occlusion of pixels within the image by shadows or parallax effects (Noh and Howat, 2015).
The present paper evaluates satellite remote sensing for the mapping of physical habitat characteristics of ice-free regions in the Antarctic Sør Rondane Mountains, in order to aid the identification and characterisation of different habitats in remote and inaccessible areas.In particular, the generation of DSM from high resolution stereoscopic imagery and derivation of Land Surface Temperature (LST) from Landsat 8/TIRS is examined.The DSM elevation data are compared against in situ recorded GPS positions and tracks, the 8 m resolution Reference Elevation Map of Antarctica (REMA) released in mid 2018 (Howat et al., 2019), and spaceborne LIDAR data.Satellite derived LST is compared against temperature loggers installed in the 2017-2018 and 2018-2019 field seasons.

Study sites
The Sør Rondane Mountains (SRM) are located about 200 km inland in Dronning Maud Land, East Antarctica, and consist of approximately 900 km 2 of nunataks with a large range of terrestrial habitats differing in bedrock and soil composition, terrain and microclimate conditions.A number of study sites (see Fig. 1 and Table 1) were selected in the SRM, based on snow cover, habitat characteristics, and accessibility from the Princess Elisabeth Antarctica (PEA) station, which operates during the austral summer season from approximately the beginning of November to the end of February.Field work was performed during January-February in three successive years, 2018, 2019, and 2020.Fig. 2 shows example photographs of the Yuboku and Pingvinane sites, with respectively large and smaller gravel fields and varying snow cover.

Satellite imagery and DSM creation
Surface elevation data are required to model the geometry of the study site, in particular the surface slope and aspect, i.e. the orientation with respect to the sun.A DSM is constructed from tristereo acquisitions from the Pléiades satellite constellation.Pléiades consists of two satellites (A/B) with on board a high resolution imaging sensor, with 4 bands (Blue, Green, Red, NIR) at 2.8 m spatial resolution, and a panchromatic channel at 0.7 m resolution.Data are resampled by the satellite operator to 2 m for the multispectral channels and to 0.5 m for the panchromatic channel.Acquisitions were made during the austral summer field work in Jan-Feb 2018, Jan-Feb 2019 and Jan-Feb 2020 (see Table 1).For each image pair in the tristereo acquisition, DSM at 1, 2, 5, 10 and 30 m spatial resolution were computed from the panchromatic channel using the Surface Extraction with TIN-based Search-space Minimization (SETSM) software (Noh and Howat, 2015, 2017, 2018).SETSM was designed especially with good performance over glaciated and mountainous regions in mind, respectively with challenging low-contrast surfaces and steep slopes and shadows.SETSM uses a Triangular Irregular Network (TIN) in object space with a pyramid of different resolutions to minimise the search space.The mean average and standard deviation of the three DSM from a tristereo acquisition were computed as the final DSM for that year.The final DSM elevation data from the last acquisition year are compared against in situ GPS positions taken at the individual study site locations, GPS tracks recorded during hikes, the 8 m Reference Elevation Model of Antarctica (REMA, Howat et al., 2019), and ICESat surface altimetry data (Schutz et al., 2005).Slope and aspect were computed from the average DSM using gdal version 2.4.1.

GPS location data
Ad hoc position and elevation measurements were performed using handheld GPS devices.Although the vertical accuracy of standard GPS is limited (about 3x the horizontal accuracy), handheld GPS devices rather than e.g.differential GPS devices were used for logging in situ elevation, due to time constraints, difficult accessibility of the sites and equipment Fig. 1.The northern part of Sør Rondane Mountains as imaged by Sentinel-2B on 2020-01-28 with study sites (cyan polygons) and sample locations (red dots) highlighted.The location of the Belgian Princess Elisabeth Antarctica (PEA) base is shown in orange.The Tvetaggen study site (north of Austkampane) is further east (at approximately 71.67 ∘ S and 25.14 ∘ E, but was not revisited in later campaigns and no data loggers were placed.Note that "Dry Valley" is a colloquial name for the rock fields between Vikinghøgda and Widerøefjellet.(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Table 1
Pléiades tristereo acquisitions made per study area for the DSM construction, with the start time provided in UTC.Sun zenith (θ s ) and azimuth (φ s ) angles are given, as well as the three viewing geometries (view zenith θ v and azimuth φ v angles).The time between first and third image in each acquisition set is about 40 s.availability.GPS data are still commonly used for recording site location and elevation, and GPS transects are included as a data source in certain elevation models (Liu et al., 1999).The measurements are included in the present study for assessing the comparability of the GPS and DSM data, and to evaluate the possibility to retrieve relative elevation transects using low-cost GPS loggers.Due to the high latitude the number of GPS satellites visible at one moment is typically very high, which should improve the provided location accuracy.GPS locations were recorded in February 2020 with a number of handheld Garmin GPS (62nd/64st/Montana 680) loggers in continuous mode for transects made at different study sites: Dry Valley (2020-02-05, 5 loggers, Fig. 3), Teltet (2020-02-06, 1 logger), Perlebandet (2020-02-07, 2 loggers) and Pingvinane (North, 2020-02-08, and South, 2020-02-09, 4 loggers).A final transect was made on 2020-02-10 at Petrellnuten (1 logger), with parts of the transect on the glacier next to the cliff wall.On a number of these transects several loggers were used in parallel to assess the vertical and horizontal errors.Logged transects were resampled to regular 30 s time grids, and when multiple loggers were active the average position on this regular time grid was computed.Additionally, per sampling location where the temperature data loggers were placed in 2018 and 2019, a fixed point was measured using the Garmin 62nd (n = 41).

Satellite LIDAR
Surface elevation data were also obtained from the Geoscience Laser Altimeter System (GLAS) on board ICESat (Schutz et al., 2005).ICESat was launched in 2003, and decommissioned after just over seven years in 2010.GLAS was a space-based LIDAR that pulses at 40 Hz and illuminates spots on the Earth's surface approximately 65 m in diameter, spaced about 172 m along track.For each of those pulses, the echo is received and a time of flight computed in order to retrieve the distance to the satellite and hence the target position.After issues with the lasers at the start of the mission, a reduced measurement program was adopted in order to track inter-annual variability of the ice sheet, with the GLAS operating three periods of 33 days each year.The GLAH12 (ice) and GLAH14 (land) products in the GLAS/ICESat L2 Global Land Surface Altimetry Dataset (version 34, Zwally et al., 2011) were obtained from the National Snow and Ice Data Center (NSIDC).Datasets covering the study site areas in Table 2 were ordered from NSIDC (accessed 2020-05-15), with elevation tracks extracted for each of the DSMs

Table 2
Boundaries (South-West to North-East) of the 12 × 12 km subscenes made per study area (plotted in Fig. 1) for the matchups with in situ temperature data.The last two columns show the number of temperature loggers retrieved in each site.(Fig. 4).For each LIDAR observation, the DSM elevation was extracted from the pixel containing the LIDAR pulse footprint, giving matchups for the Utsteinen, Dry Valley, Pingvinane, and Yuboku sites.

Surface temperature 2.3.1. Satellite imagery
Collection 1 Level 1, top-of-atmosphere terrain-corrected and orthorectified daytime Landsat 8 imagery combining data from both the Operational Land Imager (OLI) and the Thermal InfraRed Sensor (TIRS) were obtained from Google Earth Engine (Gorelick et al., 2017).OLI is a 9 band imager with 8 bands in the visible to shortwave infrared (VSWIR) at 30 m resolution and 1 panchromatic band at 15 m resolution, and TIRS has two bands (B10-B11) in the thermal infrared at around 11 and 12 μm (Irons et al., 2012).Imagery was subsetted to 12 × 12 km subscenes centred on the regions of interest where in situ data loggers were installed (Table 2).OLI data were processed using the DSF atmospheric correction (Vanhellemont, 2019) using a fixed path reflectance in the subscene in order to quality control the imagery based on the derived aerosol optical thickness (τ a ) and to derive the target emissivity based on multispectral VSWIR surface reflectance, ρ s (Vanhellemont, 2020b).

Atmospheric correction
Top-of-Atmosphere (TOA) radiance (L t ) or equivalent at-sensor Brightness Temperature (BT) is a measurement of the combined surface-emitted radiance (L s ) and atmospheric up (L u ) and down (L d ) radiances transmitted through the atmosphere (τ): To convert at-sensor BT to surface temperature, the atmospheric effects need to be corrected, for which the target emissivity (ε), atmospheric up-and downwelling radiance (L u , L d ), and atmospheric transmittance (τ) need to be determined.L u , L d , and τ are derived from radiative transfer simulations using the open source LibRadtran (version 2.0.2) radiative transfer code (Mayer and Kylling, 2005) as detailed in (Vanhellemont, 2020a, b).Two sources were used for the input atmospheric profiles of temperature and relative humidity, namely (1) a representative atmospheric profile determined from in situ radiosonde profiles (see further) and (2) the ERA5 model.Profiles were interpolated to the approximate location of PEA (71.95 ∘ S 23.34 ∘ E) at 11:00 UTC to approximate the balloon launching time (noon CET) and position.Outputs from LibRadtran were convolved to the relative spectral responses of band 10 (10 μm) and 11 (12 μm) on Landsat 8/TIRS.
Emissivity was calculated from the DSF output L8/OLI ρ s using the neural net from Vanhellemont (2020b) based on the ECOSTRESS library (Meerdink et al., 2019).

Atmospheric profiles
Atmospheric profiles were obtained from daily (2018)(2019) or every other day (all other austral summers between the 2013-2014 and 2019-2020 seasons) radiosonde launches from PEA (at around 1300 m elevation, typical pressure ± 830 hPa) around 11 UTC.Individual profiles were filtered to the upward part of the profile, and profiles were considered complete when the 99 and 1 percentiles of pressure were respectively <700 hPa and >200 hPa.Profiles with a 99 percentile pressure >850 hPa were also excluded as anomalously high at PEA. Temperature and humidity measurements of complete profiles were interpolated between 900 and 0 hPa pressure at regular intervals of hPa.Complete profiles (n = 243) were averaged across each field season and across the full dataset.Pressure values were filtered if the number of samples in the average was less than 10, and pressure was subset every 50 hPa for input into LibRadtran.The complete radiosonde profiles for the 2018-2019 season and seasonal averages are shown in Fig. 5.

In situ temperature
Small (16 mm) iButton Hygrochron temperature and relative humidity loggers (Maxim Integrated, USA) were installed in the SRM during the 2018 and 2019 field campaigns, configured to record temperature (T, ± 1 ∘ C and relative humidity (RH, %) every 3 h.The manufacturer quotes the iButton uncertainty as ± 1 ∘ C, but in certain conditions results better than 0.3-0.5 K are obtained (Brewin et al., 2019(Brewin et al., , 2021)).These accuracies are of a similar magnitude to the noise-equivalent temperature difference of the TIRS instrument (0.8 K at 240 K and 0.4 K at 300 K, Irons et al., 2012), but smaller than the expected performance for general land temperature matchup performance (unbiased RMSD around 3 K, Vanhellemont, 2020b).As both satellite and in situ measurements have a relatively large uncertainty, statistics are described as differences between both datasets, and the Reduced Major Axis regression lines are used to assess their relative performance.The logging settings were chosenwith memory and battery limitations in mindto allow for continuous operation for a one-year deployment.Loggers were installed on the upper surface of the bedrock, covered with a small rock to reduce the risk of loss, and to improve recoverability.After approximately one year the loggers were retrieved during the following field season.43 and 39 loggers were retrieved with data spanning respectively 2018-2019 and 2019-2020.An example of in situ measured time series of temperature is given in Fig. 6 for two loggers placed only metres apart in the Dry Valley site.Two snow fences were installed in 2018 to artificially increase snow cover after the loggers were put in place.Note the large temperature difference between the exposed and snow-covered logger, and the large temporal variability during the polar day of the exposed logger.The snow cover also leads to lower annual variability (see also results from e.g.Convey et al., 2018), with a temperature range observed here between approximately − and − 10 ∘ C for the snow-covered logger, and a range between approximately − 40 and +15 ∘ C for the exposed logger.
Matchups between satellite and in situ measurements were automatically identified from the daytime 12 × 12 km L8 subscenes with DSF derived τ a at 550 nm ≤ 0.25, and availability of bounding three-hourly in situ measurements.Bounding in situ measurements were linearly interpolated to the satellite overpass time, and the data range between the bounding measurements was stored.The satellite data were taken as the mean average of a 3 × 3 pixel box from the 30 m OLI grid, centred on the deployment location.The standard deviation from the 3 × 3 pixels was also stored.Matchups with in situ data range <2 K in the 3 h window were classified as likely frozen or snow covered (e.g.Fig. 6), while matchups with larger in situ data ranges experienced some direct heating by the sun.

Supervised classification
To perform a supervised classification per site, the Landsat derived products are composited in time, filtering out cloudy and hazy scenes.The average, standard deviation and maximum LST derived from Landsat 8/TIRS and surface reflectance derived from Landsat 8/OLI were computed and reprojected to the 10 m DSM data extent.The maximum LST is chosen to represent a typical summer high surface temperature, and is computed for each pixel in the Landsat composites.The Normalised Difference Snow Index (NDSI) was computed based on the Landsat 8 mission average surface reflectance (ρ s ) in the Green (560 nm) and SWIR1 (1609 nm) bands: NDSI = ρ s 560 − ρ s 1609 ρ s 560 + ρ s 1609 . (2) Each pixel was classified according to the NDSI and DSM derived variables.Four classes of snow/ice cover were derived from the NDSI: NDSI <0: no ice, NDSI 0-0.4 sparse ice cover, NDSI 0.4-0.8moderate ice cover and NDSI≥0.8fully ice covered.The DSM derived slope was separated into four classes: < 5 ∘ , 5-15 ∘ , 15-25 ∘ , ≥ 25 ∘ , the aspect in five classes: the four cardinal directions (North, South, East, West ± 45 ∘ ) and slope < 5 ∘ (i.e.low to no slope, and hence no real orientation), and the elevation in four classes: < 1100m, 1100-1300 m, 1300-1500 m, ≥ 1500 m.Pixels were extracted for different class combinations, and the  mean average maximum LST was computed from the mission composites.

Digital Surface Model
An example of the Digital Surface Model (DSM) constructed from the Pléiades tristereo imagery of the Pingvinane Nunataks is presented in Fig. 7.The retrieved DSM elevation data are compared with spaceborne LIDAR data, the 8 m Reference Elevation Model of Antarctica (REMA) made available in mid 2018, and in situ measured GPS positions and tracks (Fig. 8).The Pléiades DSM elevations agree well with the ICESat/ LIDAR observations made a decade or more before, with a mean bias of − 6.5 m MAD, i.e. the Pléiades DSM giving lower values, and a scatter of about 8.8 m RMSD.Similar values are found for the REMA with − 6.2 m MAD and 9.2 m RMSD.REMA was generated using the same processing software and use of similar optical imagery (GeoEye-1 and WorldView-1/2/3), acquired during the austral summers of 2009-2017, with the bulk of the imagery from 2016 to 2017.REMA was vertically registered to data obtained from CryoSat-2 radar and ICESat GLAS during the 2-D campaign at the end of 2008 (Howat et al., 2019), and validated to have a vertical error < 1 m compared to airborne data.
In order to evaluate the GPS accuracy at the latitudes of the SRM, GPS data from 5 handheld devices were recorded during a climb up the Northern part of the Dry Valley site on 2020-02-05 (Fig. 3).Measured positions were interpolated to a regular grid in time (every 30 s), and the mean and standard deviation across the 5 loggers was computed.The standard deviation was generally within 2 × 10 − 5 degrees for latitude, and within 7 × 10 − 5 degrees for longitude, corresponding to horizontal errors less than approximately 2.3 m at 72.1 ∘ south.Elevation ranged from 1600 to 1950 m, with elevation errors between 1 and 4 m (5-95 percentiles) and median error 2.8 m.Different performance was found for the ascent and descent, perhaps due to the difference in speed: the ascent took about 1 h, and the descent 30 min.Even though GPS elevation data have a large vertical errortypically about 3x the horizontal position errorthe differences between the DSM and GPS are reasonably small (MAD 16.4 m and RMSD 24.8 m, or about 1-2% of the target elevation).Although the vertical accuracy of GPS is known to be relatively poor, GPS point and transect measurements are common sources of elevation for hard to reach areas, and are shown here to be rather consistent across several loggers, presumably as a result of the low moisture content of the atmosphere (humidity delays the GPS signal travel time and hence the overall accuracy, e.g.Bevis et al., 1992) and high number of concurrently observable GPS satellites over the Antarctic.A significant part of the errors encountered in the matchup between DSM and GPS are coming from the Petrellnuten elevation track, where part of the track was made on the glacier in the wind scoop in close proximity to the cliff protruding >100 m above the ice.For these pixels the differences are very large due to mixed pixels, positioning errors, and occluded visibility from the satellite.Such effects can also be observed in the satellite dataset for other very steep features (e.g.wind scoop ice walls).Excluding the Petrellnuten dataset reduces the errors to MAD of 12.3 m and RMSD of 15.2 m.
The accuracy of the DSM presented here is limited because only the satellite rational polynomial coefficients were used, and no bias  Pléiades DSM generated here can be shifted to reduce the biases if higher elevation precision is needed, e.g. using better ground control data, or one of the reference datasets.This kind of bias correction was done for REMA (Howat et al., 2019), but has not been performed here for the current application.

Atmospheric profile impacts
Profiles of atmospheric temperature from the PEA radiosondes were fairly stable (Fig. 5), and even though some variability was present in the atmospheric profiles of relative humidity, the variability in the seasonal averages was low.The modelled atmospheric parameters based on the annual average profiles vary little, and have quite low radiances and high transmittances due to the low heat capacity of the atmosphere: L d 0.069-0.111W/m 2 sr − 1 μm, L u 0.039-0.062W/m 2 sr − 1 μm, τ 0.985-0.990,for B10 excluding the 2014-2015 season.The 2014-2015 season had on average a quite warm atmosphere (see bump at around 750 hPa in Fig. 5), leading to lower atmospheric transmittance (τ 0.802) and higher atmospheric radiances (L d 2.511 W/m 2 sr − 1 μm and L u 1.541 W/m 2 sr − 1 μm).The average over all seasons gave for B10: L d 0.216 W/ m 2 sr − 1 μm, L u 0.119 W/m 2 sr − 1 μm, and τ 0.976.Slightly larger variability is found across the 243 individual complete profiles, 5-95 percentiles for B10: L d 0.045-0.211W/m 2 sr − 1 μm, L u 0.026-0.118W/ m 2 sr − 1 μm, and τ 0.975-0.993.Although relative humidity profiles above the tropopause may not be reliable due to potential earlier condensation or ice build-up on the sensor, the impact of these errors is very low.Excluding the upper part of the profile, with pressure <200 hPa, leads to differences of <1% in the estimated atmospheric parameters, and mean differences of <0.001K in the final LST product.Differences between the radiosonde and ERA5 profiles are relatively small, for B10, on average around 0.1 K (with 0.15 K standard deviation) for L t of 5 W/m 2 sr − 1 μm and − 0.1 K (with 0.10 K standard deviation) for L t of 6 W/m 2 sr − 1 μm (Fig. 9).Estimated LST differences between the PEA profile and the average PEA profile show a similar mean average difference, but a narrower histogram (i.e.smaller standard deviation).These differences are well below the Noise Equivalent Temperature Difference of the TIRS specification, 0.8 K at 240 K and 0.4 K at 300 K (Irons et al., 2012), and hence the ERA5 model is a suitable source of atmospheric profiles for processing Landsat 8 LST data.The error in LST will not depend on the atmospheric profiles, but will be largely determined by small scale variability (snow, ice, and rock cover) around the data logger positions, with some smaller impacts of the sensor itself and the estimated per pixel emissivity factor.

Surface temperature matchups
For the 2018-2019 and 2019-2020 seasons, respectively 500 and 401 matchups were identified with the in situ iButton data from 5 different sites.A summary scatter plot is given in Fig. 10 for both seasons (901 matchups), the per scene derived emissivity, and the ERA5 atmospheric profiles, giving an overall Mean Difference (MD) of − 2.5 K and an unbiased Root-Mean Squared Difference (unb-RMSD) of 6.3 K.The unusually warm conditions in December 2019 are causing the cluster of outlier data in the bottom right of the plot (around 0 ∘ C in situ and − 25 ∘ C satellite temperatures), where there is a clear mismatch between the in situ point measurement and the satellite pixel.By including data only from before December 2019, 818 matchups remain, with a MD of − 1.8 K and an unb-RMSD of 4.6 K (not shown separately).Matchups were classified as either frozen/snow covered or exposed based on the in situ temperature range between the logged temperatures bounding the satellite overpass.The absolute values of the MD are overall smaller for the matchups classified as frozen, i.e. − 1.8 K and − 1.3 K, than for the matchups exposed to solar heating, i.e. − 4.0 K and − 2.7 K, respectively including and excluding December 2019 data.This means that a large part of the bias can be attributed to temperature changes between the two bounding measurements as a result of solar heating, and that higher frequency in situ logging may reduce the bias.The bias for the stable frozen points is close to what was observed in Belgian sites with the same processing system (Vanhellemont, 2020b), with higher scatter (unb-RMSD).The scatter for the present results is rather consistent between the two classes of points, i.e. 5.8K and 7.0K for the data including December 2019, and 4.6 K for both classes excluding December 2019.These results also indicate that the scatter in the results is likely influenced by the sub pixel scale (100 m) variability around the logger deployment location, which can vary on the centimetre scale due to rock and ice coverage, shadow and direct sun exposure (Convey et al., 2018).There could be some additional differences caused by the difference between the emitted radiance derived temperature and the temperature measured by contact temperature probes, and the accuracy of the temperature probesquoted as ± 1 ∘ C by the manufacturer for the − 30 to +70 ∘ C range.
Performance was found to be different per site, consistent over both campaigns (Table 3).For example, for sites with flat gravel and rock beds extending over multiple 100 × 100 m TIRS pixels (Yuboku and Dry Valley) a positive MD is found, i.e. satellite retrievals are warmer than in situ measurements.This is caused by the satellite averaging the top of the gravel and rocks that are illuminated and heated by the sun, and the data loggers measuring under the top gravel layer, generally fixed with a pebble to improve recoverability after one season in the field.At Yuboku for example, during the austral summer, shallow meltwater lakes may be present due to the preferential heating of the rocks and gravel, and the shape and orientation of the valley, even though air temperatures remain well below freezing.Similar meltwater ponds and liquid water on rock surfaces and in damp soils have been observed even further south, at Ellsworth Land 75-76 ∘ S (Convey and McInnes, 2005) and the Dufek Massif at 82 ∘ S (Hodgson et al., 2010).For the other sites, negative MD are found, i.e. satellite retrievals are colder than in situ measurements, with generally − 4-− 5 K MD for the narrow rock outcrops (Utsteinen, Perlebandet) and about − 1.5 K for the nunatak slopes (Pingvinane).For these sites, rock and gravel make up only a small portion of the TIRS pixel, and hence the satellite data contains large parts of the surrounding, significantly colder, ice and snow fields.Overview photos of the Yuboku and Pingvinane sites are shown in Fig. 2.
For the PEA radiosonde profiles data, 473 matchups (52%) had an available profile taken within one week, but this reduces the data to the field season timeframe, i.e. − 20 ∘ C and warmer.For the data with radiosonde profiles, the overall MD was − 2.8 K with unb-RMSD of 7.8 K (note that this includes the December 2019 data).For the ERA5 profiles, the same subset gave very close results: MD of − 2.9 K and unb-RMSD of 7.9 K.An alternative processing using the closest season average PEA radiosonde profile for the matchups with no profile taken within one week allows for processing of the full 901 matchup dataset, and gives near identical errors compared to the ERA5 dataset: MD of − 2.5 K and unb-RMSD of 6.3 K.For large scale automation over multiple areas, the use of ERA5 profiles will hence provide adequate results.The use of ERA5 will provide results for sites or periods without available balloon data -about half of the matchups in the case of PEA.

Classification results
The remote sensing datasets could help explain e.g.differences observed between soil and microbial community samples taken at different sites.This is not further explored in the current manuscript, but a classification method is demonstrated here that may be useful for future analyses.Fig. 11 shows the classification results for the Pingvinane nunataks based on the multispectral OLI data (NDSI) and the DSM derived characteristics.LST data can be extracted for any combination of these classes, for example, the ice-free regions have much higher temperatures and variability compared to the glaciers, and it is found that East and North facing slopes generally have higher temperatures than the West and South facing slopes (examples not shown separately).
The Landsat 8 mission composites clearly show the expected temperature differences between snow and ice covered areas and the icefree areas, with an example for Yuboku Valley shown in Fig. 12. Maximum temperatures differ by a few degrees between the white snow packed glacier and the snow-free blue ice situated close to the mountains, with the blue ice being about 2 ∘ C warmer.In Yuboku Valley, the ice-free areas reach maximum LST of 10-15 ∘ C and higher, especially for the East and North facing mountain slopes and the central valleywhere the Yuboku lakes are located.The West facing slopes are generally shaded during the Landsat overpass time (between 6:42 and 7:02 UTC) and hence show much lower temperatures compared to the East and North facing slopes, with differences in the maximum temperature reaching up to 10-20 ∘ C.
Analysis of data from temperature loggers placed at Pingvinane North (on the West face) and Pingvinane South (on the North face) indicate that there is a clear difference in maximum temperature between the North and West facing sites at the time of Landsat overpass, but that there is not necessarily a difference in maximum daily temperature at the sites (Fig. 13).The temperature peak on the West face lags the North face by about 3 h, and in fact these measurements indicate that at Pingvinane, the West facing sites experience higher temperatures compared to the North facing sites.It is hence clear that the morning overpass of Landsat 8 leads to strong biases in the LST data, especially since the area receives 24 h of sunlight of varying intensity during the polar day.These results indicate that a more complete modelling of the solar irradiance and surface heating based on the DSM is needed for a more detailed physical habitat analysis of the ice-free regions.The Landsat LST could still provide valuable validation data for this analysis.

Conclusions
Our results demonstrate the potential of satellite remote sensing for characterising surface temperature and site elevation and orientation of Antarctic ice-free mountain areas at metre and decametre scale.Elevation was retrieved from Pléiades and SETSM with an RMSD of 5 m, and temperature was retrieved from Landsat 8/TIRS with a MD of <2 K and RMSD of <5 K. Attention should be paid to the local overpass time of Landsat and the impact on the satellite derived LST.Due to the limited range of overpass times of Landsat there are significant biases in the mission average temperatures -even though the matchup results show biases of only a few degrees.This is no different from other regions on earth, but since there generally is 24 h of sunlight during summer in the Antarctic, the temperatures recorded at the Landsat overpass are not at all representative of the daily maxima reached for the different mountain faces.Since the temperatures encountered in the Antarctic ice-free regions are often around the freezing point of water, the difference of a few K of the LST compared to the in situ measurements may have significant impacts if the LST data is used to model biological processes, and has to be interpreted with care.However, L8/TIRS does provide an unprecedented high resolution surface temperature dataset for the Antarctic ice-free areas.
Strong evidence is presented that the temperature matchup performance is influenced by the in situ data logging frequency, and also the site specifics, i.e. the spatial coverage of rock/gravel compared to snow/ ice within a 100 × 100 m TIRS pixel.Sites with consistent gravel/rock composition over several TIRS pixels showed different responses compared to sites where rock made up only a small fraction of the pixel.When the rock composition dominates, the satellite gives higher temperature due to their measurement of the top of the rocks heated by the sun, whereas the in situ loggers measure under the upper layer of gravel.Conversely, when the ice fraction dominates, the satellite data give lower temperatures compared to the in situ loggers installed on the (warmer) rock outcrops.Our results indicate that higher spatial resolution thermal data would be required to fully resolve the heterogeneous temperature encountered in the Antarctic ice-free regions (e.g.Convey et al., 2018).While the TIRS-2 onboard Landsat 9 (planned for launch in September 2021) will be nearly identical to the L8/TIRS, the planned Land Surface Temperature Monitoring (LSTM) mission of the European Union's Copernicus programme may provide higher spatial resolution TIR imagery.The LSTM mission requirement document (Koetz et al., 2019) lists 30 and 50 m spatial resolution (or approximately 0.5 and 1 ha) as goal and threshold requirements respectively.While these future missionsand physical limitations of thermal sensor designwill not be able to resolve separate habitats on a single nunatak slope, they will continue the L8/TIRS temperature data record.A large temporal variability of temperature due to solar heating is found within a 3 h window for non-ice covered data loggers, and a linear interpolation to the overpass time may underestimate the temperature.Higher frequency logging, or a more appropriate interpolation model e.g.incorporating non-linear heating effects, may further reduce the bias between satellite and in situ measurements.Future work should include a thorough assessment of the iButton performance at the low temperatures encountered in the SRM (− 40-20 ∘ C), and impacts of snow and ice cover.Potential further applications of the high resolution DSM include the computation of the seasonal light availability from direct sun and diffuse sky irradiance, which may influence local temperature evolution.From these solar irradiance simulations, the heating potential could be estimated for different mountain sides and rock types for the full daily and seasonal cycles.The Landsat LST product could then be used for validation of these modelled temperatures, and perhaps rock types could be derived from the OLI mission composite surface reflectances and a collection of reference spectra (e.g.ECOSTRESS).If higher resolution or higher accuracy DSMs are needed, e.g. to study centimetre scale terrain effects, other techniques such as drone overflights or ground-based laser scanning are needed to complement the metre and decametre scale

Fig. 2 .
Fig. 2. Photos showing Yuboku overview (top left) and lake (top right) and Pingvinane overview (bottom left) and iButton in gravel before covering with pebble (bottom right).

Fig. 4 .
Fig. 4. Tracks of the ICESat/GLAS data through the SRM and the 12 × 12 km study sites.Of these the Perlebandet site was the only one where tracks did not intersect the Pléiades DSM.Colour intensity relates to the number of tracks in the GLAS archive.(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 5 .
Fig. 5. Temperature and relative humidity profiles as obtained from noon radiosonde launches at PEA.Individual complete profiles from the 2018 and 2019 field season (top), and summer average profiles across the different field seasons (bottom).

Fig. 6 .
Fig. 6. iButton temperature time series for one of two snow fences installed for the MICROBIAN project at the Dry Valley study site.Two loggers were placed <5 m apart, one on exposed gravel (blue) and one was covered under the snow pack created by the snow fence (orange).(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 10 .
Fig. 10.Left, matchups between in situ recorded temperature and Landsat 8/TIRS B10 derived temperatures, using the ERA5 atmospheric profiles, and per pixel OLI derived emissivity.Horizontal error bars show the data range between the two three-hourly in situ measurements that bounded the satellite overpass, and vertical error bars are the spatial variability (3 × 3 pixel standard deviation) as recorded by satellite.The right hand side plot shows the results separated for likely frozen or snow covered (blue), and non-frozen (orange) loggers.(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 11 .
Fig. 11.Classification of the Pingvinane site according to Normalised Difference Snow Index, and DSM derived Elevation, Slope and Aspect.

Fig. 13 .
Fig. 13.Maximum temperatures recorded for iButtons placed at Pingvinane on a North facing (5 loggers) and West facing (2 loggers) slope for ten days in November 2019.The red vertical lines denote the typical Landsat 8 overpass time (around 7 UTC).Average slopes are 10 ∘ and 14 ∘ respectively for the North and West facing slopes, both with average elevation around 1415 m. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Table 3
Results of temperature matchups per season and study site.December 2019 data was excluded.