Analysis of Lava from the Cumbre Vieja Volcano Using Remote Sensing Data from DESIS and Sentinel-2

Abstract

On 19th September 2021, a protracted eruption of the Cumbre Vieja Volcano on the Canary Island of La Palma commenced and continued for a duration of 12 weeks. Lava flows starting from the rift zone at the mid-western flank of Cumbre Vieja advanced toward the western coast of the island. The eruption was monitored by different remote sensing satellites, including the Copernicus Sentinel missions and DESIS. The Sentinel-2 Copernicus satellites acquired multispectral data from 15th September onward. On September 30th, and with a difference of ∼2 h with respect to Sentinel-2 A, the DESIS hyperspectral sensor also acquired data from the volcano and then again on 15th October 2021. Typically, mid-infrared (around 3.8 $\mathsf{\mu }$m) data are used for the thermal analysis of active lava flows. However, neither Sentinel-2 nor DESIS possesses mid-infrared bands and the Sentinel-2 high-wavelengths bands (∼2 $\mathsf{\mu }$m) have some limitations. Nevertheless, the hyperspectral character of DESIS enables the analysis of active erupting volcanoes in near-infrared wavelengths. The results of this analysis find fluid lava temperatures of about 1100–1200 K but there are problems associated with the high-temperature lava spectral emissivity.

1. Introduction

During large, rapidly evolving volcanic eruptions, when it is difficult to gain in situ information due to safety reasons, satellite remote sensing plays a key role in supporting the monitoring activities. Thereby, analysis of thermal remote sensing data enables the investigation of a variety of thermal volcanogenic emitting phenomena, such as lava lakes [1], active lava flows [2,3], and lava domes [4]. The thermal remote sensing of volcanoes has been applied since the early 1980s, beginning with the Advanced Very High Resolution Radiometer (AVHRR) onboard the National Oceanic and Atmospheric Administration (NOAA) satellites and the Landsat Thematic Mapper (TM) series [5,6].

In particular, the joint investigation of mid (MIR, around 3.8 $\mathsf{\mu }$m) and thermal infrared (around 10 $\mathsf{\mu }$m) via the dual-band method [7] was used to estimate lava temperature and area [8,9]. To retrieve more accurate information of the lava flow, the dual-band approach has been extended further (e.g., assuming three thermal components) [10]. Nowadays, automated processors such as MODVOLC [11] and MIROVA (Middle InfraRed Observation of Volcanic Activity) [12], both based on the Moderate Resolution Imaging Spectrometer (MODIS), enable the detection and monitoring of thermal volcanic anomalies worldwide. MIROVA also provides information about the volcanic radiative power (VRP), directly derived from MIR-band data using the approach of [13], which is valid for hot objects with temperatures above 600 K and provides information on the VRP with an uncertainty of ±30%.

Satellite-data-based analysis of the temperature regime of an active lava flow (as well as the shortcut of Wooster’s MIR-band approach [13]) enables estimates of the time-averaged discharge rate (TADR), the lava volume flux averaged over a given time period [14]. The TADR is very important as it is a major controlling factor of the lengths of a lava flow and therefore also a critical parameter for predictive models of lava flow aerial development [15]. Eruptions of Stromboli Volcano in 2014 [16,17] and 2019 [18] show significant overestimation in the satellite-based estimations when compared to the topographic change detection. But, depending on the type of eruption, and especially on the thickness of the developing crust, thermal data alone have high uncertainties regarding the absolute lava volume, but provide valuable information on short-term fluctuations in the lava effusion rate [4].

More detailed analysis of effusive volcanic eruptions is only possible with higher-spatial-resolution sensors such as Sentinel-2 or Landsat-8/9 (with 10 to 30 m GSD) [19]. These sensors do not have an MIR channel but provide information in the short-wave infrared (SWIR) and near-infrared (NIR) spectral regions. This article investigates the estimation of lava flow temperatures from high-spatial-resolution SWIR and NIR data, analyzing observations from multispectral Sentinel-2 and hyperspectral DESIS (DLR Earth Sensing Imaging Spectrometer).

By focusing on optical remote sensing acquisitions in the visible to near-infrared (VNIR) to SWIR wavelengths range (Section 2), new findings can be reported regarding the emissivity of lava at magmatic temperatures. During this 12-weeks long eruption, a total volume of 217.4 ± 6.6 × 10${}^6$ m${}^3$ was emplaced, comprising 36.5 ± 0.3 × 10${}^6$ m${}^3$ from the volcanic edifice and 177.6 ± 5.8 × 10${}^6$ m${}^3$ from subaerial lava flows (including fallout deposits on the lava flows, but excluding the submerged portion of the lava deltas) [20], along with a tephra blanket of approximately 20 × 10${}^6$ m${}^3$, accounting for only 7–16% of the overall erupted volume [21,22].

As a case study, we investigate the large 2021 lava flow event at Cumbre Vieja Volcano, La Palma, first analyzing Sentinel-2 multispectral data (Section 3.1) and later adding more spectral information using hyperspectral data (Section 3.2). Discussion about the analysis results and conclusions will be addressed in Section 4 and Section 5.

2. Optical Remote Sensing Data

The analysis of the temperature radiated by ground objects is typically analyzed from the SWIR (1 $\mathsf{\mu }$m < $\lambda$ < 3 $\mathsf{\mu }$m) to the thermal infrared (8–14 $\mathsf{\mu }$m), depending on where the maximum of the blackbody emission is located. According to Wien’s displacement law, Earth surfaces with ambient temperatures around 300 K emit a maximum radiation close to the 10 $\mathsf{\mu }$m wavelength. With increasing surface temperature, the peak emission shifts to smaller wavelengths. As an example, for T = 800 K, the peak emitted radiance is around 3.7 $\mathsf{\mu }$m, and, for T = 1100 K, at 2.7 $\mathsf{\mu }$m (compare Figure 1). The 3.7 $\mathsf{\mu }$m spectral region is located in an atmospheric window region and the radiation can be recorded by Earth-observing remote sensing satellite sensors with channels in this region, such as MODIS, VIIRS (Visible Infrared Imaging Radiometer Suite), FireBird, and Sentinel-3 SLSTR (Sea and Land Surface Temperature Radiometer).

For objects with even higher temperatures (1000–1200 K), the emitted radiation can also be recorded in shorter-wavelength regions (at approximately 2.2 $\mathsf{\mu }$m) (see Landsat-8-9 [23,24,25], Sentinel-2 [26], EnMAP (Environmental Mapping and Analysis Program ) [27], PRISMA (PRecursore IperSpettrale della Missione Applicativa) [28], EMIT (Earth Surface Mineral Dust Source Investigation) [29], among others). High temperatures can also be measured with channels around 1.6 $\mathsf{\mu }$m or around 0.7–1.0 $\mathsf{\mu }$m, even if these channels are not located at the peak emission wavelength. In the NIR/SWIR wavelength range, a complication exists for daytime data, where the reflected and scattered solar radiation has to be subtracted from the total signal when calculating the radiation temperature. If the surface emissivity is known or estimated, then the surface kinetic temperature can be easily retrieved.

Figure 1 presents examples of radiation spectra for three Earth temperatures (ambient 300 K, high at 800 K and 1100 K) together with the reflected solar radiation for three surface albedos ($\rho =0.01,0.1,1.0)$. Here, the influence of the atmosphere is neglected. In the 10 $\mathsf{\mu }$m region, the solar reflected radiance can be neglected, but for wavelengths <4 $\mathsf{\mu }$m, both components have to be taken into account. In particular, for wavelengths smaller than 2 $\mathsf{\mu }$m and ambient temperatures $(\sim$300 K), the solar component dominates.

However, for optical sensors taking data in the NIR/SWIR wavelength domain (<2 $\mathsf{\mu }$m), the radiation emitted by high-temperature events (700–1500 K) can be larger than the solar component. These sensors can then detect the radiation emitted by high-temperature ground objects. In particular, the Sentinel-2 bands B8A (865 nm), B11 (1600 nm), and B12 (2200 nm) can be used for this purpose (see Section 3.1). For the DESIS hyperspectral sensor, more bands are available but only in the VNIR range (see Section 3.2).

For this investigation, we employ the IDL software package ATCOR (ATmospheric CORrection) [30] and its Python successor PACO (Python-based Atmospheric COrrection) [31]. Both software packages characterize and correct the atmosphere for optical remote sensing sensors. For their atmospheric correction, they use previously simulated atmospheric conditions, stored as look-up tables (LUTs). These radiative transfer functions will be used for atmospheric variables in this analysis.

Energy Balance Equation

By the conservation of energy, all the incident radiation in a surface is partly transmitted (TR), partly reflected (R), and partly absorbed (A). Thermal equilibrium is reached when a body emits the same energy that it absorbs (A = E). For opaque objects, TR = 0. This leads to the final relation between radiation fractions: R + E = 1. Optical remote sensing is a passive technique that measures the reflected solar irradiance as the source of energy and/or the thermally emitted radiance. The emitted part allows for the calculation of the brightness temperature, and even the surface temperature if the emissivity is known. For daytime data acquisitions, objects illuminated by the sun typically have surface temperatures < 350 K, and thermal equilibrium can be assumed. The emissivity $\epsilon$ is then related to the surface reflectance $\rho$ via $\epsilon \left(\lambda \right)=1-\rho \left(\lambda \right)$.

However, during a volcanic eruption, the active hot lava emits an additional radiation source independent of the solar warming, and the thermal equilibrium with respect to the solar irradiance is not fulfilled. Therefore, reflectance ($\rho$) and the emissivity ($\epsilon$) of the body are considered independent, and, for high temperatures, the emissivity also depends on the temperature [32,33]. As a summary, for daytime data, the satellite signal in the region 800–2500 nm consists of reflected solar radiation, solar path radiance, and emitted radiation by high-temperature sources.

In this scenario, we neglect the solar path radiation in the 700–2500 nm spectrum because this component is very small (typically < 1% for $\lambda >$ 850 nm) compared to the reflected components (direct and diffuse). Thus, the total top-of-atmosphere (TOA) radiation L, received at the sensor, can be written as:

$$L=L_r+L_e$$

(1)

where $L_r$ is the reflected direct solar radiation (Equation (2)) and $L_e$ (Equation (3)) is the emitted radiation.

For the solar component of the energy equation, the diffuse solar radiation on the ground can be neglected in the SWIR (1000–2500 nm) since it is typically only 10% (at 1000 nm), 5% (at 1600 nm), and 3% (at 2200 nm) of the direct solar component.

The reflected solar radiation can then be expressed as:

$$L_r=\rho E_{0}cos{\theta }_s{\tau }_s\tau /\pi$$

(2)

where $\rho$ is the surface reflectance, $E_0$ is the extraterrestrial solar irradiance, ${\theta }_s$ is the local solar zenith angle, ${\tau }_s$ is the atmospheric transmittance of the sun-to-ground path, and $\tau$ is the transmittance of ground-to sensor path.

Neglecting the thermal ambient emission from the body due to the solar component, the emitted spectral radiance can be expressed as:

$$L_e(\lambda ,\epsilon ,T)=\epsilon \left(\lambda \right)L_{bb}(\lambda ,T)$$

(3)

where the object emissivity dependency of the temperature for such high temperatures [33] is included in $\epsilon \left(\lambda \right)$ and the emitted blackbody spectral radiance $L_{bb}$ depends on wavelength $\lambda$ and temperature T according to Planck’s equation,

$$L_{bb}(\lambda ,T)=L_e/\epsilon =\frac{2h{c}^2}{{\lambda }^5}\frac{1}{exp(hc/\lambda k_{B}T)-1}$$

(4)

where h is Planck’s constant, c is the speed of light, and $k_B$ is the Boltzmann constant.

The solution of the energy balance equation involves the determination of three variables: the surface reflectance $\rho$ (or scene background $L_r$, Equation (2)), the body emissivity ($\epsilon$), and the temperature (T). The analysis approach and results for multi- and hyperspectral data are discussed next.

3. Analysis and Results

The combination of Equations (1)–(3) approximates the remotely sensed TOA radiance acquired by different sensors. Additionally, the acquisition channels of these sensors have a certain bandwidth and the resulting radiance equations must be convolved with the spectral channel response function $R\left(\lambda \right)$ of a sensor integrated over the bandpass.

Furthermore, since the measurements are not at ground level, the atmospheric transmittance $\tau$ (corresponding to the sensor view angle) and ${\tau }_s$ (for the solar geometry) will be taken into account. The ground elevation per remote sensing pixel (via the digital elevation model, DEM) and the visibility (in km) of the atmospheric conditions are also considered, extracting the sun-to-ground and the ground-to-sensor transmittance functions from radiative transfer (RT) simulations stored in look-up tables (LUTs).

The recorded radiance data per band can be compared with the model using two different approaches:

• 1. Extract the temperature per band. With only one available NIR channel in the case of Sentinel-2 (the SWIR channels were saturated), the emissivity of the hot body at the band central wavelength is assumed and the reflectance is extracted from the adjacent ground. Then, the remaining variable, the temperature of the hot body, can be determined. This approach will be discussed in Section 3.1.
• 2. Model the radiance through a least squares fitting of the unknown variables ($\rho$, $\epsilon$, and T) together. This approach can only be carried out when more than three spectral channels are available. This is the case for hyperspectral sensors (Section 3.2).

In order to compare the results from both methodologies, the analysis will be applied to Sentinel-2 and DESIS data during their overpass over the lava flow on 30th September 2021.

3.1. Sentinel-2

The Sentinel-2 A scene (GRI tile 28 RBS) was acquired on 30th September at 13:00 UTC (Coordinated Universal Time). Figure 2 (left) contains an overview of La Palma Island containing the red lava flow area (west coast of the island), and the right part shows the cosine of the local solar zenith angle, derived from the solar zenith angle (34.2${}^{\circ }$), solar azimuth angle (154.5${}^{\circ }$), and digital elevation model (DEM) [34].

Figure 3 is a zoom of the west coast in Figure 2 containing the lava flow. Figure 4 contains the zoomed view of the lava region for different band combinations.

Sentinel-2 spectral data consist of thirteen bands at different wavelengths and three different spatial resolutions (10 m, 20 m, 60 m) [35]. In order to make a complete spectral analysis of the Sentinel-2 results, we have merged the spectra of all the Sentinel-2 bands into a single spectral cube. The spectral bands of the 10 m and 60 m data cubes have been interpolated spatially and merged to a final 20 m resolution TOA data cube with 13 spectral bands.

Among the 13 bands, Sentinel-2 has 2 bands in the SWIR region (named B11 and B12) and 1 in the NIR region (B8A). Figure 4 (top) shows an RGB color composite with band B12. A comparison with the B11 composite (middle) clearly shows the strong blooming effects for band B12, indicating detector saturation in the hot lava region. This fact is not so obvious for band B11 but, nevertheless, saturation also occurs for most lava pixels, as demonstrated when analyzing their spectra in detail (see Figure 5).

Analyzing the 20 m TOA reflectance spectra, DN of the scaled radiance, and TOA radiance from lava location (x,y) = (911,1561) using a 3 × 3 pixel average on Figure 5, we see very high values in TOA reflectance ($\rho >1.2$) in bands B11 (1600 nm) and B12 (2200 nm). The original ESA encoding of the TOA reflectance scene is 16 bits/pixel using a scale factor of 10,000 [26]. So, the value $\rho =1.4$ for band B11 (1600 nm) in Figure 5 (top left) is coded as DN = 14,000 in the TOA reflectance image cube, far below the maximum 16 bit value. In the ATCOR and PACO models, the TOA reflectance is converted into scaled TOA radiance, labeled DN in Figure 5 (top right), with channel-dependent scale factors. The DN and radiance values are also detailed in

Table 1. DN and TOA radiance (mW cm${}^{-2}$ sr${}^{-1}$$\mathsf{\mu }$m${}^{-1}$).

Channel	DN	DN (bg)	L (TOA)	L (TOA, emitted)
B8A	8430	1462	8.430	6.968
B11	44,086	3380	8.817	8.142
B12	60,048	15,012	3.002	2.252

. In addition, the typical background radiance (“DN (bg)”) from a 3 × 3 pixel neighborhood close to the lava flow is subtracted from the total radiance to obtain the radiance emitted from the hot, fluid (i.e., not crusted) lava. The corresponding physical radiance is given in Figure 5 (bottom), clearly showing that not only is band B12 saturated but also band B11.

Therefore, in the Sentinel-2 scene (La Palma, 30th September 2021), the high radiance emitted by the lava flow area in the SWIR bands is above the maximum radiance of the dynamic range ($L_{max}$) for those bands (69.78 and 24.6 W m${}^{-2}$ sr${}^{-1}$$\mathsf{\mu }$m${}^{-1}$ for B11 and B12, respectively) [36] (see Table 1), so we cannot evaluate temperatures here.

The bottom image in Figure 4 contains a band combination with channel B8A (865 nm). There are notably fewer “red” lava pixels in the lava area (the other red pixels outside the lave field indicate vegetation) and when the corresponding radiance and temperatures are evaluated (Figure 5 and Table 1), it is clear that no saturation is encountered in band B8A (865 nm).

Table 1 shows the digital number (DN) of the scaled TOA radiance, the background DN(bg), and the corresponding TOA radiances for bands B8A, B11, and B12. The emitted TOA radiance is calculated with Equation (5). It accounts for the subtracted radiance $L_{bg}$:

$$L(TOA,emitted,i)=L(TOA,i)-c1\left(i\right)DN(bg,i)$$

(5)

Then, we assume $\epsilon =0.97$ for all channels [33] and $\tau =0.85,0.90,0.90$ (B8A, B11, B12, respectively) to calculate the surface radiance (Equation (6)) from the emitted TOA radiance and obtain surface temperatures for the fluid lava of 1118 K, 764 K, and 573 K for bands B8A, B11, and B12, respectively (

Table 2. Lava temperatures: T1 = TOA, T2 = TOA-background, T3 = surface.

Channel	T1 (K)	T2 (K)	T3 (K)
B8A	1117	1103	1118
B11	760	755	764
B12	580	566	573

). In the case of $\epsilon =0.90$, we obtain a slightly higher B8A lava surface temperature of 1123 K for the same L(TOA, emitted).

$$L_{bb}\left(T\right)=L(TOA,emitted)/\left(\tau \epsilon \right)$$

(6)

This clearly indicates that not only band 12 (B12) but also B11 is above the maximum radiance of the dynamic range. Unfortunately, this is not coded as the maximum of the 16 bits/pixel encoding in the TOA reflectance cube. Therefore, only the B8A TOA temperature of 1117 K is valid (T1 in Table 2). Next, the background radiance is subtracted from the total radiance (T2 in Table 2). The last step is the calculation of surface temperatures (T3 in Table 2).

Comparing the T1 and T2 temperatures demonstrates the negligible influence of the background radiance. A comparison of T1 (TOA) and T3 (surface) shows that the atmospheric transmittance (and, to a lesser degree, the assumed emissivity) compensate for the background radiance influence, and the temperatures T1 and T3 are almost the same. The B11 and B12 values are not reliable because of the detector saturation.

The temperature results for the non-saturated NIR band B8A are compatible with the hyperspectral results detailed in Section 3.2. Since the Sentinel-2 data have only one non-saturated band available for a quantitative evaluation, assumptions about the emissivity and surface reflectance were necessary for the temperature evaluation. The next section calculates results for a hyperspectral sensor (DESIS), where more spectral information is acquired, avoiding the assumptions about the emissivity and reflectance.

3.2. DESIS

The above problem concerning the saturation of the SWIR bands does not exist in the VNIR spectrum of DESIS data. Since many bands are available, the three unknown parameters (temperature, emissivity, surface reflectance) can in principle be retrieved with a least squares fit algorithm with sufficient degrees of freedom (number of bands). This is one of the advantages of hyperspectral sensors as detailed below.

DESIS [37] is a VNIR sensor mounted on board the International Space Station (ISS) and it has been operational since October 2019. DESIS acquires images of the Earth at user request with a swath of about 30 km, 30 m spatial resolution, and 235 bands with a full width at half maximum (FWHM) of 3.5 nm in the spectral range from 400 to 1000 nm.

The DESIS data used in this analysis are the L1C product, the ortho-rectified remote sensing measurements as TOA radiance (in mW cm${}^{-2}$ sr${}^{-1}$$\mathsf{\mu }$m${}^{-1}$).

The DESIS overpass of the La Palma Island took place approximately 2 h after the overpass of Sentinel-2 on 30th September 2021, and again on 15th October.

Table 3. DESIS data acquisition conditions of 30th September and 15th October.

Date	Time	datatakeID	$\mathbf{\Theta }$ (${}^{\mathbf{\circ }}$)	Off-Nadir (${}^{\mathbf{\circ }}$)
30th September	14:12	638761160	35.6	19.3
15th October	17:01	644308832	69.4	5.4

summarizes the overpass conditions of the datasets, shown in Figure 6 as RGB quicklooks.

Only bands in the NIR wavelengths (>700 nm) will be used in the model fit.

3.2.1. Hot Fluid Lava Mask

The first step in the analysis consists of the selection of the set of pixels that are most likely part of the hot, fluid lava flow. The rest of the pixels (>90% of the image) will not produce convergent results in our analysis, so their consideration is neglected for computing power purposes. Therefore, only a reduced number of pixels belonging to a pre-selected “lava mask” (focusing on the hot fluid parts of the lava flow) will be used in this analysis.

By adapting the normalized hotspot index (NHI) approach proposed by [38,39] to the DESIS NIR bands (NHI${}_{NIR}$), this lava mask is calculated using the ground apparent reflectance (${\rho }_{app}$) for the following condition:

$$NH{I}_{NIR}=\frac{{\rho }_{app}\left(995\mathrm{nm}\right)-{\rho }_{app}\left(850\mathrm{nm}\right)}{{\rho }_{app}\left(995\mathrm{nm}\right)+{\rho }_{app}\left(850\mathrm{nm}\right)}$$

(7)

The spectral increase due to the thermal influence is selected by NHI${}_{NIR}$> 0.1 and water fluctuations due to noise are excluded by ${\rho }_{app}\left(995\mathrm{nm}\right)>0.1$.

3.2.2. Model Fit to Hyperspectral Data

Through a least squares fitting of the DESIS hyperspectral L1C data, we obtain per pixel the T, emissivity ($ϵ$), and reflectance ($\rho$). All three parameters are assumed as constant within the wavelength range between 700.0 and 999.5 nm, where the model fit residuals converge.

As for Sentinel-2, the atmospheric transmittance (sun-to-ground, ${\tau }_s$, and ground-to-sensor, $\tau$) is extracted from MODTRAN simulations [40], convolved with the DESIS sensor spectral response functions for the sun zenith angle and off-nadir angle of the observations, the corresponding pixel elevation above sea level (DEM), and an atmospheric visibility of 23 km.

Even inside the pre-selected lava mask, some pixels do not have a large thermal radiance component to give good convergence results. These pixels are excluded from the results. For this analysis, we only considered those pixels where the chi-square of the least square fit chi-square obtained is above the 99% confidence level for a chi-square distribution with 3 free parameters and 39 degrees of freedom (d.o.f), assuming an uncertainty in the DESIS L1C radiance of ∼4% [41].

Figure 7 shows the fit results for one of the lava pixels of DESIS. The DESIS pixel (x = 413, y = 747) corresponds geographically to the position analyzed in the Sentinel-2 (Section 3.1) scene of the same day. For this pixel, the obtained fit has an emissivity $ϵ$ = 0.12 ± 0.01, T = 1291 ± 11 K, and a reflectance $\rho$ = 15.9 ± 0.5%.

Once the procedure is verified with the test pixel, it is applied to the full set of pixels in the lava mask.

The statistics over the resulting active lava flow are summarized in

Table 4. Mean and standard deviation ($\sigma$) of the active lava flow parameters (temperature T, emissivity $ϵ$, and reflectance $\rho$) (weighted mean and standard deviation) from DESIS acquisitions.

Date	T [K]	RMSET [K]	$\mathit{ϵ}$	RMSEϵ	$\mathit{\rho }$ [%/100]	RMSE${}_{\mathit{\rho }}$ [%/100]
30th September	1159 ± 2	111	0.09 ± 0.001	0.25	0.17 ± 0.00	0.02
15th October	1237 ± 2	184	0.040 ± 0.001	0.246	0.260 ± 0.001	0.1

for the DESIS scenes of 30th September and 15th October.

The values in the table represent the weighted mean and standard deviation of the total active (fluid) lava flow area. These results assume that all pixels belonging to the active lava flow share the same properties. Therefore, the fit to each pixel is a measurement of the same element and the outliers are due to less accurate fit results. Table 4 gives the mean and its error averaged by the inverse of the pixel fit error for each parameter.

However, the lava flows are not necessarily homogeneous and they can present a gradient distribution along their paths, decreasing the temperature with the distance to the eruption center. Since the reflectance and the emissivity might also change along the lava flow, the population distribution of the measurements of any of the three fitted parameter cannot be considered Gaussian. Therefore, Table 4 gives the distribution mean of the temperature, reflectance, and emissivity and their RMSE. The distribution of the active, fluid lava has a mean temperature of 1134 K, emissivity of 0.32, and reflectance of 17%. The minimum/maximum values (min, max) of the resulting distribution of the three lava flow parameters are: T [K] = (897, 1398), $ϵ$ = (0.01, 0.97), $\rho$ [%/100] = (0.10, 0.24).

If we select the pixel of the DESIS scene with the lowest ${\chi }^2$ value, (x, y) = (389, 742), we obtain a temperature of 1129 ± 6 K, an emissivity of 0.32 ± 0.02, and a reflectance of 16.7 ± 0.2%.

Figure 8 shows the map of the fitting results per pixel for the emissivity (top), temperature (middle), and reflectance (bottom) for the DESIS scene of 30th September. Pixels with emissivity values of 1 are actually fit results that have reached the boundary limit of the emissivity parameter. For the scene of 15th October, only five pixels obtained good fit parameters, so results cannot be displayed as a map. The low number of pixels might be due to the amount of clouds around the lava flow.

4. Discussion

Due to the limited number of bands from multispectral sensors as well as the decrease in the LFWC (linear full-well capacity) in the SWIR bands for some optical remote sensing instruments, an insufficient number of bands (or even just one band, B8A, in the case of Sentinel-2) are available for the determination of the three parameters of the energy radiance balance in the NIR. The reason, of course, is that most high-spatial-resolution satellite sensors are not designed to measure high-temperature events.

Even in the case of hyperspectral sensors like DESIS, which has many more bands, the thermal signal at high wavelengths is quite marginal (see Figure 7).

Nevertheless, the values extracted for the temperature of the active lava (∼1200 K) seem to be compatible with other independent measurements (∼1090–1214 ${}^{\circ }$C) [42]; however, this compatibility is only indicative since the temperature depends on the location (the RMSE of our lava flow results already has a spread of ∼100–200 K), the time of measurements, and whether the lava is fluid or already partially crusted at the surface.

The constant emissivity of ∼0.1–0.3 between 700 and 1000 nm is completely outside the published emissivities for these hot objects [33,43], where the wavelength-dependent values are typically between 0.9 and 1. The resulting values from the fit in the NIR wavelengths do not fulfill the relation $ϵ=1-\rho$. Further night observations could help to validate the emissivity results obtained since, there, the solar irradiance component is discarded as input in the energy balance equation, and only the emissivity and the black-body emission contribute to the total radiance received at the sensor.

The limited available DESIS wavelength range (<1000 nm) is far from the wavelengths (∼2–3 $\mathsf{\mu }$m) at which the peak of the thermal emission is expected. Other models [43], accounting for sub-pixel sources and more than one temperature, do not give convergent results. These models, where the thermal emission is limited to a fraction of the pixel area, could increase the emissivity by reducing the emitting area from 1 to 0.3–0.5.

Independent measurements of the typical emissivity of this type of lava in the NIR wavelength could also be conducted and then added as input to the fitting model, helping to further constrain the temperature results.

The addition of SWIR bands into the analysis (multi- and hyperspectral) is often limited by the full-well capacity (FWC) or the maximum radiance of the dynamic range of the sensor at those wavelengths. The limited amount of charge that can be stored within an individual pixel without the pixel becoming saturated (definition of the FWC or LFWC) is a design feature and most Earth remote sensing satellite instruments are not designed for measuring high-temperature events. This feature becomes more important in high-spatial-resolution sensors (e.g., 20, 30 m spatial resolution), where the detector dynamic range of one or more pixels can be completely saturated by the high thermal component of a lava flow. This is in contrast to sensors with >100 m resolution, where the high thermal component usually only covers a certain fraction of the dynamic range of a detector pixel.

The analysis of the lava flow from L1B or L1C remote sensing products presented here takes, at the moment, 2–3 min. This time does not include any parallelization optimization, since the amount of pixels inside the lava flow is rather small (197 out of 2.4 × 10${}^6$ pixels in the image). To be used for volcano crisis management, parallelization processes per pixel can and should be implemented. The final processing time should be added to the mean L1B/L1C products access time for the remote sensing mission.

5. Conclusions

The 2021 Cumbre Vieja eruption, La Palma, was analyzed by the high-spatial-resolution remote sensing missions Sentinel-2 (multispectral sensor) and DESIS (hyperspectral). On 30th September 2021, they observed an active lava flow with a time difference of 2 h.

After a limited analysis due to the DESIS wavelength range available for the measurements ($700<\lambda <1000$ nm), the temperature results along the active lava flow yield ∼1000–1300 K. These values are of the same order of measurements with Sentinel-2 band B8A (865 nm). However, the obtained constant emissivities between 700 and 1000 nm are much lower (∼0.1–0.3) than the typical values at short-wave infrared wavelengths (with typical emissivities ∼0.9–1) and they clearly do not fulfill the relation $ϵ=1-\rho$. Further in situ validation studies and analysis results on night observations could be useful for determining the possible limitations of the determination of the emissivity using NIR wavelengths.

Using the normalized hotspot index (NHI) adopted to DESIS data (Equation (7)), we detected in the DESIS image acquired on 30th September fluid lava at the top third (upstream) part of the lava flow, with the hottest temperatures (up to 1398 K) being observed within the center line of the lava flow and lower temperatures at its rim. The vents area was affected too much by ash clouds to detect fluid lava there, while the lower two-thirds of the lava flow were characterized by lower temperatures at the surface due to crust development, which isolated the hotter fluid lava flowing in tubes beneath it. Larger cracks in this crusted surface are shown by “islands” of fluid lava further downstream within the center of the lava flow. For 27th September 2021, i.e., 3 days prior to our first DESIS acquisition, a strong increase in the lava effusion rate is described [4]. Based on seismic data indicating the activation of a deeper magma source (reactivation of the 12–15 km cluster) on 27th September, the authors assume a change in the lava composition from tephrite to basanite, i.e., to less viscose lava, as was observed in previous eruptions at La Palma (e.g., in 1949 [44] and 1971 [45]). At least for the upstream part of the lava flow, where we detected fluid lava in the DESIS image of 30th September, the less viscose basanite lava played a dominant role. However, in order to derive information about the lava viscosity directly from the DESIS data themselves, a time series of data would have been more helpful. Analysis of high-spatial-resolution hyperspectral data with SWIR wavelengths is recommended for further studies provided that the SWIR channels do not saturate.

But, the analysis of active lava in the NIR spectrum is a challenge because of missing emissivity measurements at these high temperatures. Therefore, additional emissivity measurements of this kind of lava at magmatic temperatures at NIR wavelengths will help to validate these results and reduce their uncertainty, and can be used to set realistic boundaries for fit parameter models. Finally, wavelength-dependent lava emissivity models should be included in future temperature evaluation studies of multispectral and hyperspectral sensors whenever the corresponding data are available.