Inverting multispectral thermal time-series images of volcanic eruptions for lava emplacement models

Satellite thermal images have been used to observe volcanic activity for several decades (e.g. Gawarecki et al. 1965; Hanel et al. 1979; Rothery et al. 1988; Glaze et al. 1989; Oppenheimer & Francis 1997) and have found particular use in estimating different ‘rates’ of emplacement of volcanic material on planetary surfaces, such as volumetric effusion rate, time-averaged discharge rate (TADR) and rates of areal emplacement (Harris et al. 2007a). Such estimates have been used to infer the behaviour of the magmatic systems generating volcanic activity (Harris et al. 2000, 2003; Ripepe et al. 2005), parameterize lava flow emplacement models (Wright et al. 2008; Vicari et al. 2009), find relationships between rates of emplacement and other geophysical parameters (e.g. Coppola et al. 2009), and place constraints on extraterrestrial volcanic activity (Carr 1986; Davies 1996).

There are two main ‘classes’ of models in the literature relating thermal radiance and rate of lava emplacement, which we will refer to as the instantaneous radiance models and the temperature distribution models. The instantaneous radiance models, developed over a series of papers (Harris et al. 1997a, b; Wright et al. 2001; Harris & Baloga 2009), assume that the TADR is proportional to the active area of the lava flow, based on the analytical model and empirical relationships of Pieri & Baloga (1986), with a constant of proportionality that can be derived from the model (despite some contradictions pointed out by Dragoni & Tallarico 2009) or estimated empirically. The active flow area is found by inverting radiances for a two-component subpixel thermal model (Dozier 1981).

The temperature-distribution approach is based on the observation that as a lava-flow advances it exposes fresh hot material that cools in a relatively simple and predictable way. This cooling trajectory combined with an area coverage rate defines a surface temperature distribution, which in turn defines an emission spectrum (Carr 1986). The emplacement history is effectively ‘written in’ to the lava surface temperature distribution, and the emission spectrum at a point in time is therefore a function of the area emplacement rate up to that point. This basic conceptual model is supported by field observations for a number of different lava-flow types. For instance, for channelized lava flows, the downflow temperature profile typically exhibits an exponential decrease consistent with regular cooling of crust (Oppenheimer 1991; Pinkerton et al. 2002), and converting distance from vent to crust age using the flow speed gives cooling curves that can be modelled using simple finite-difference models (Ball et al. 2008), while pahoehoe lobes at varying stages of emplacement exhibit temperature distributions that ‘decay’ in a regular fashion, as evidenced by the histograms presented by Wright & Flynn (2003). In both cases, there is a relationship between the temperature of a surface, its age and, therefore, its emplacement history.

One approach to utilizing these observations to recover emplacement rates is to attempt to recover directly the temperature distribution itself. If thermal images of the lava flows are of sufficient resolution to capture the spatial variation in temperature, a direct ‘thermal chronometry’ (Harris et al. 2007b) approach can be taken by relating imaged temperature variation to age variation (e.g. Davies 2003; Harris et al. 2007b). Otherwise, emission spectra have to be inverted to give the subpixel temperature distribution, a procedure known as the ‘inverse blackbody problem’ (Sun & Jaggard 1987). Most terrestrial approaches use the methods of Dozier (1981) and Oppenheimer (1993) to model the pixel as hot incandescent cracks, cooler crust and cold background; however, it has been demonstrated that around six thermal components are necessary to adequately describe the temperature distributions of lava flows found from high-resolution thermal camera surveys (Wright & Flynn 2003). The limited number of bands provided by most sensors and the constraints placed by atmospheric windows limit the number of unknowns that can be recovered, and, typically, some temperatures and fractional areas have to be assumed a priori. Even in the absence of an atmosphere, and with approximately 540 bands between 4 and 55 µm, Pearl & Sinton (1982) were only able to model thermal spectra of Io as a function of three separate Plank functions. In addition to these practical difficulties, the inverse black-body problem is fundamentally ‘ill posed’, limiting the extent with which the temperature distribution can be resolved (Lakhtakia & Lakhtakia 1986).

Given the frequent subpixel nature of the spatial variation in temperature, the complex shape of volcanic temperature distributions, and the practical and theoretical difficulties associated with the inverse black-body problem, attempts have been made to model temperature distributions and emission spectra and time series as a function of a few emplacement parameters and the lava thermophysical properties (Carr 1986; Davies 1996). Some approaches have assumed a constant effusion rate of infinite duration, which defines a time-invariant equilibrium thermal-emission spectrum, while others assume constant effusion rate but assume it starts and stops in time, defining an emission spectrum that varies with time. The former only allows for inversion of emission spectra for area coverage rate; however, the latter also allows inversion of a time series of observations at a given wavelength, in addition to the inversion of spectra. Furthermore, some investigations have attempted to take into account cracking and subduction of the surface crust area into the lava flow, and overriding of the lava surface at the front of the flow, by including a constant fraction of any area element radiating at a predicted core temperature (Davies 1996), or by defining a survival time for a lava-flow surface area before it is recycled into the flow (Carr 1986; Keszthelyi & McEwen 1997). The overturn of lava lakes or seas has been modelled by varying resurfacing distributions with stationary resurfacing fronts, migrating resurfacing fronts or foundering (Matson et al. 2006).

Here we attempt to generalize these models to interpret the emplacement and removal of lava surface at any temperature, while still assuming that all surface elements follow the same cooling trajectory. This gives the radiance in wavelength and time as a function of a Net Area Emplacement (NAE) surface that varies with time and surface emplacement temperature convolved in the time direction and integrated in the temperature direction with a radiance cooling kernel. In other words, this gives a universal model of the radiant output of any style of volcanic activity where the surface of erupted material cools in a similar manner, which will be more applicable to lava lakes and gentle effusive eruptions than vigorous fountain-fed lava flows. The NAE surface contains useful information about the nature of the processes adding and removing hot material in a sensor’s field of view and the overall rate of expansion of the anomalously hot region, and so is a very desirable retrievable for thermal remote sensing of volcanic activity. However, the overlapping nature of the radiance observation kernel functions and typically limited number of bands make the equation underdetermined and ill posed, fundamentally limiting the resolution that can be attained, and resulting in a wide family of NAE surfaces that provide an equivalent fit to the observed radiances. We present a simple regularized inversion scheme and test it on combinations of typical instrument wavebands and plausible NAE models to explore the parameters that could be reasonably retrieved given an assumed cooling curve.

We first present the derivation of the model and show how it can be discretized to give a simple linear matrix equation. We then show how the nature of the integration kernels and limited sampling in wavelength result in an underdetermined and ill-posed problem, and how we can attempt to stabilize the solution with regularization. We end by describing the plausible artificial NAE models we use to test the inversion technique.

Equation (7) is linear and so is amenable to the techniques of linear inversion, which allow investigation of its mathematical properties and inversion. This requires discretizing the equation in three variables: wavelength, time and temperature or temperature-equivalent age, so the equation can be approximated as a discrete matrix equation. In the following we discretize assuming that we have thermal anomaly radiance values recorded at regular time intervals over the same geographical region and at the same wavelengths, a situation that only applies for geostationary orbiters and stationary terrestrial radiometers. More sophisticated discretizations could be used to describe radiances acquired irregularly in space, time and wavelength by a variety of low Earth orbit instruments.

Given an assumed kernel matrix K, we aim to invert the radiance vector R for the NAE vector A.

The kernel matrix K in equation (20) is primarily a function of the surface area element cooling curve, T(t). The cooling curve can be calculated using models of cooling lava in different geometrical configurations relevant to particular modes of emplacement, or it can be taken from physical measurements (e.g. time series of thermal camera images, radiometer measurements or thermocouple measurements) – here we adopt the modelling approach. Two geometries can account for a wide range of volcanic phenomena: lava cooling in a column can be used to simulate emplacement of successive areas of an advancing lava flow, while lava cooling as a sphere can simulate the cooling of scoriae. Both scenarios can be modelled using simple finite-difference methods based on the Taylor expansion of the one-dimensional (1D) heat equation (Özisik 2002), which have been used extensively in various forms to model the cooling of lava (e.g. Kent & Pinkerton 1994; Keszthelyi & Denlinger 1996; Ball et al. 2008). In the columnar scenario, the lava is modelled as a vertical column of elements all at an initial ‘eruption’ temperature, T_e, and heat is lost from the surface by radiance and convection, and transmitted through the column by conduction; in the spherical scenario, however, the lava is modelled as a series of concentric shells, the outer of which loses heat by radiation and convection, with heat transmitted between shells by conduction. For more details and a full derivation, see Özisik (2002).

In this study we use an adaptive time step and mesh size to speed up computation, increasing step and/or mesh size when the temperature gradient of the surface element with respect to space or time drops below some threshold. The mesh size and time step are updated subject to the constraint that they remain within the stability field of the finite-difference scheme. The model produces a series of (time, surface temperature) pairs that can be interpolated to give surface temperature as a continuous function of time. Example runs of the two models are shown in Figure 2.

The radiance from the surface element is simply the Planck function of the temperature and the relevant wavelength. However, in this study we are attempting to model thermally anomalous radiance (i.e. radiance in excess of that expected in the absence of an eruption), so we subtract the expected radiance from the ambient temperature from that predicted by the model in order to compensate for the missing radiance due to the lava surface obscuring the background surface during emplacement.

We have discretized equation (7) in terms of wavelength bands and samples every Δt in the time direction and at points t_0,i in the temperature-equivalent-age direction in order to get the matrix equation (equation 20), which describes the forward problem, or the calculation of a radiance vector R from a kernel matrix K and a model vector A. The tractability of inverting R for A is dependent upon the mathematical properties of the matrix K, and we can analyse the nature of K in order to understand the limitations imposed on the inversion by the physics of the model using the standard toolkit of linear inversion, which is summarized in a number of textbooks (e.g. Twomey 1996; Hansen 1998; Gubbins 2004; Aster et al. 2011). This is principally a function of the kernel vectors (Twomey 1996), or rows of the matrix K, which, taken with the dot product of the model vector, give each radiance value, one per row: in other words, each row of K gives the dependency of the corresponding radiance observation on the NAE surface. We can take the kernel vector for a given radiance observation and ‘unpack’ it to show it as a function over the NAE surface, as shown by the examples in Figure 3.

Figure 3 shows kernels for observations 900 s apart in time and at wavelengths in the shortwave infrared (SWIR), mid-infrared (MIR) and thermal infrared (TIR), a typical sampling interval and set of wavebands for geostationary imagers (e.g. SEVIRI (Spinning Enhanced Visible Infra-Red Imager): Aminou 2002). Four kernels are shown for observation 50 in a time series of 100 at four different wavelengths, over a range of typical volcanic temperatures, for both columnar and spherical lava geometries. The kernels vary with wavelength and observation number: for a given observation time, all show significant overlap with varying wavelength, and, for the same wavelength, all show some overlap in the observation direction, particularly for the longer wavelengths. This spread of the kernels results in variations in the NAE being averaged out over the area of the kernel, resulting in small changes in observed radiance, making radiance a weak function of detail in the NAE. The overlap of the kernels results in them being strongly correlated, resulting in small variations in differences between radiances for large variations in the NAE. This means that, for the inverse problem, small changes in radiance can be mapped to large changes in the NAE surface, rendering any attempt at inversion unstable and sensitive to noise. It is this correlation between kernels that makes the inversion ill posed and our kernel matrix K ill conditioned (Twomey 1996). The fact that the variation in the sensitivity of the kernel to NAE with increasing wavelength manifests by increasing overlap rather than offset (as is the case with varying observation number) limits the utility of increasing the number of bands per observation (Twomey 1996).

The problem is usually underdetermined in that the matrix K is rank deficient, which may be due to: (i) fewer observations than unknowns since multispectral images typically collect only a few bands that can register volcanic radiance, giving a few radiances per value of t, while we may desire to invert for several t₀ values per t; and (ii) bands being so close together in wavelength that some effectively do not contribute extra information, resulting in rows of K that are linear functions of each other. This results in the inverse having a non-trivial ‘null space’, or a number of vectors in the NAE vector space that are mapped to the zero vector in the radiance vector space, such that there are certain ways in which the NAE vector can vary that give no change in the radiance vector: for instance, certain patterns of area emplacement at higher temperatures and removal at lower temperatures may cancel each other out in radiance terms. In other words, there are NAE vectors that cannot be differentiated on the basis of the patterns of thermal radiance that they produce.

We use equation (24) to model an NAE on a grid that consists of 100 observations at intervals of 900 s in the time direction (which is the acquisition interval for SEVIRI, and a typical order of magnitude for geostationary satellites) and at temperature intervals of 30K in the T/t₀ direction, an example of which is shown in Figure 4. For the cooling curve we use a columnar model with an eruption temperature of 1375K, ambient temperature of 300K, thermal diffusivity of 7×10⁻⁷ m² s⁻¹, thermal conductivity of 1.5 W m⁻¹ K⁻¹, emissivity of 0.95 and convection coefficient of 60 for a 2 m vertical column of lava. We model two types of NAE surfaces, which we refer to as ‘simple’ and ‘complex’. In both scenarios, we model material being emplaced with positive NAE values across a range of temperatures to simulate a priori the effects of view angle and diverse geomorphological processes on initial exposure temperature within the sensor IFOV (instantaneous field of view), and with a skew in the time direction to simulate the waxing and waning effusion rates frequently seen in fissural eruptions. In the ‘complex’ scenarios, we add a similar negative NAE feature at lower temperatures to simulate the removal of material. In both cases, the two NAE features are moved in the NAE plane with temperature, as varying the models in time or in magnitude should have a comparatively minor effect. The NAE used here has units of (pixel fraction) K⁻¹ s⁻¹.

Radiances are then calculated from these NAE surfaces using equation (20), and Gaussian noise with a standard deviation of 0.05 times the standard deviation of the radiance signal is added. We invert for combinations of radiances from bands at 1.6, 3.9 and 10.8 µm, which we refer to as SWIR (shortwave infrared), MIR (mid-infrared) and TIR (thermal infrared) to give potential combinations for common satellites (e.g. SEVIRI), and also a combination that contains all SWIR, MIR and TIR bands plus two in the NIR (near-infrared) (0.8 and 1.0 µm) to test the effect of adding extra shorter wavelengths. An example of the inversion process is shown in Figure 4. The effectiveness of the inversion is examined visually for some examples, and in terms of the retrieval of specific NAE features (temperature of maximum NAE, total NAE, total positive NAE and total negative NAE) as a function of this shift in temperature.

Example inversions for two different ‘simple’ models are shown in Figures 5 and 6. The top row of panels shows the test model NAE surface and the TAE curve, which is just the NAE integrated over the t₀/T direction to give the net rate of change of area as a function of time: multiplying the latter by a mean material thickness would give a volumetric effusion rate. The pairs of surface and line plots below show the recovered NAE surface and TAE curve for a given combination of wavebands. In both Figures 5 and 6, all waveband combinations smear out the NAE surface to some degree, with single-band solutions smearing it across the full range of temperatures, and increasing numbers of bands better localizing the positive feature, giving a smaller support and higher peak NAE value. In all cases, the NAE is well constrained in the time direction. Shorter wavelengths, alone or in combination, tend to over- or underestimate the TAE curve the most, while longer wavelengths tend to reproduce it most faithfully. Interestingly, despite smearing the NAE surface badly, TIR alone, and TIR and MIR, do a good job of reproducing the TAE curve, an observation we will come to focus on later.

We might expect the quality of the recovery to vary most with location of the positive NAE feature in the temperature direction due to the aforementioned overlapping nature of the kernels in this direction (as opposed to the offset nature in the t direction), so we ran the inversion for NAE surfaces with the mean of the feature shifted in 100K increments, and tested the recovery by finding the percentage difference in the total amount of material emplaced (integral under the whole NAE surface) between the modelled and recovered surfaces, and temperature at which the NAE is at a maximum in the modelled and recovered surfaces. The results are shown in Figure 7. For high temperatures, the shorter wavelengths overestimate the total, and for lower temperatures they underestimate the total substantially (by up to 50%); however, most other band combinations are within 5–10%. Single- and double-band combinations locate the maximum NAE value at either the top or bottom of the model domain, swapping between the two as the positive feature moves with temperature. Three or more bands are sufficient to localize the peak somewhere in the middle of the domain, but the location is frequently out by up to 100K.

An example inversion for a ‘complex’ model is shown in Figure 8. Again, the top row of panels shows the test NAE surface and TAE curve, and the panels below the results of the inversion. In this case, we have a negative feature offset at lower temperatures to represent the removal of material from view by, for example, rollover, subduction or view geometry effects, with the same dimensions but a smaller amplitude to give a positive TAE curve. Again, single-waveband inversions smear the NAE across the whole temperature range of the domain, while increasing numbers of bands progressively localize the positive and negative features. Shorter wavelengths alone and in combination over- or underestimate the total emplacement curve, while TIR alone, TIR and MIR, and three or more band combinations reproduce it well, albeit with some noise at sharp features that manifests as spurious ‘wiggles’ due to the Laplacian regularization. We again test the effect of displacing the model in the temperature direction, and test the effectiveness of: (i) the recovery of the total positive part of the NAE surface (the total area emplaced, ignoring that removed); (ii) the total negative part of the NAE surface (the total area removed, ignoring that emplaced); (iii) the total amount emplaced (integral under the whole NAE surface, both positive and negative parts); and the recovery of the temperatures at which the NAE is at (iv) a maximum and (v) a minimum. Each of these tests is shown as a panel in Figure 9.

We are particularly interested in the total positive and negative parts of the NAE surface, as isolating them tells us something about the amount of resurfacing taking place, independent of the net expansion of the hot surface. The plots show that the total positive and negative components of the NAE surface are consistently underestimated, with the best results from three or more band solutions underestimating the total positive by 20–30% and total negative by around half. The total area (integral over both positive and negative NAE features) is again poorly recovered by short wavelengths alone or in combination, and well recovered by the TIR alone and multiple wavebands. As before, the maximum and minimum NAE values tend to be pushed out to the edges of the domain, with only the five-band solution effectively localizing anything in the middle of the domain.

Finally, we performed two last inversions on complex models with positive and negative NAE features that are close in total value such that the net total emplaced area is small: this can be considered to model a situation in which there is a lot of ‘churn’ at the hot surface and most material is removed in some way. In Figure 10, the NAE has a the same support at high temperatures as the previous ‘complex’ examples, and, as before, shorter wavelengths alone and in combination smear the NAE and badly misestimate the total area emplacement curve, while the TIR more faithfully reproduces the curve. The TIR and MIR pair, and the three or more band combinations, reproduce the curve but with a substantial amount of noise, again manifest as spurious ‘wiggles’, and also localize the NAE more effectively. In Figure 11, the NAE has a broad support across the range of temperatures in the domain resulting unusually in a poor curve recovery for the TIR band alone, but a reasonable, if noisy, recovery for the TIR and MIR pair, and three or more band inversions.

The test models presented do not exhaustively sample the full space of possible NAE surfaces; however, we have explored a range of plausible scenarios from which we are able to draw broad conclusions about which aspects of the NAE can be recovered using different waveband combinations. The quality of the recovery of the NAE surface and the TAE curve are primarily a function of the number of bands and the shape of their kernels with respect to the domain upon which the NAE surface is defined.

In general, recovery of the NAE and the TAE curve tends to be better in the time direction than in the temperature direction, as expected given the offset nature of the kernels in the time direction and the overlapping nature in the temperature direction. Good resolution in the time direction is important, as post-eruption cooling of hot material can be mistaken for continuing emplacement. Increasing the number of bands increases the number of unknowns and improves recovery, with combinations with longer wavelengths performing better.

Single-band solutions do not provide enough unknowns to resolve the NAE in the temperature direction, smearing it out across the whole domain as a function of the balance between minimizing the model residual and its curvature. However, bands with a relatively ‘flat’ kernel across the temperature direction conserve NAE with smearing such that the TAE curve is accurately retrieved regardless. For instance, the TIR kernel in Figure 2 is relatively flat across the full temperature range, such that smoothing an arbitrary NAE surface subject to maintaining the same radiant output requires the addition of the same ‘amount’ of NAE at one location as is removed at another. This accounts for the good recovery of the TAE curve in Figures 5, 6, 8 and 10. In contrast, the SWIR kernels in Figure 2 vary substantially across the domain in the temperature direction, such that smoothing an arbitrary NAE surface subject to a given radiant output does not conserve NAE. For instance, the best-fit solution can replace a small quantity of NAE at high temperatures with a large quantity at low temperatures and maintain the same radiant output. In fact, at the lowest temperatures the SWIR kernel is effectively zero, and the surface is unconstrained and free to vary solely for minimizing curvature, resulting in the retrieval of positive NAE and TAE values before eruption onset in Figures 5, 8, 10 and 11.

The recovery of the curve is improved moderately if the MIR band is added, although, in the case where the curve is the small difference between two large positive and negative NAE features, this can significantly amplify noise (e.g. Fig. 10) which the TIR-only inversion is not subject to as the kernel is averaging across temperature. Increasing the number of bands improves the resolution of the NAE surface, localizing the support of the surface, and sometimes localizing the temperature of maximum and/or minimum NAE within the domain rather than at the edges, and improving the estimate of the magnitude of the positive and negative features of the NAE surface. However, even for the best inversions with three or five bands, localized NAE maxima are frequently out by 100K or more, and the magnitude of total emplacement and removal are underestimated by up to 60%

We have presented a novel generalization of the family of area emplacement models frequently used to interpret multispectral time series of thermally emitted radiance from volcanic eruptions. The equation describes the time series of radiances in terms of a convolution-integration of two functions, a Net Area Emplacement (NAE) function and a cooling curve. The NAE describes the emplacement or removal of surface area at any time and temperature or temperature-equivalent age, while the cooling curve describes the cooling of an area element with time. The generalization rests upon the crucial assumption that all surface elements follow the same cooling curve, independent of the length of exposure or initial temperature to which they are exposed in the instrument FOV.

We have presented a discretization of this equation to frame it as a simple matrix equation, shown how, owing to the shape of the kernels for each radiance observation, the inverse problem is ill posed, and shown how it can be inverted using a common regularization technique. Testing the inversion with plausible test NAE models for combinations of typical wavebands used in orbital sensors has allowed us to draw some broad conclusions, namely: (i) three or more wavebands spanning the SWIR, MIR and TIR allow recovery of the broad shape of the NAE, although quantification of features such as the balance between total emplacement and total removal are heavily biased; (ii) the inversion successfully accounts for cooling post-emplacement despite these biases; and (iii) inversion of single wavebands can recover the total area emplacement curve (the rate at which the hot area is expanding) well if the kernel is flat across the range of temperatures in the NAE domain. These conclusions all have implications for the thermal remote sensing of volcanic eruptions and suggest future work, as follows.

First, the main assumption of the model, namely that most of the material emplaced follows a similar cooling curve, needs to be tested in a range of settings. There are a number of eruption types where the dominant emplacement processes cannot all be approximated by the same kernel: for instance, the radiant signal from lava flows fed by large fire fountains would need to be modelled using both the spherical and columnar kernels shown in Figure 2. The model probably applies reasonably well to lava lakes and gentle effusive lava flows; however, real NAE surfaces need to be found to compare with our ‘plausible’ a priori parameterized models to check whether our deductions are reasonable. Such NAE surfaces could then be used as test models themselves to refine our understanding of the biases introduced by the inversion process.

Another area where the model could be improved is in the inversion itself. The inversion is heavily underdetermined and ill posed, and here we resolved this by finding a minimum curvature solution; however, alternative regularization techniques may do a better job of suppressing the noise seen in the recovered total area emplacement curves in Figures 10 and 11. In addition, it may be that extra constraints can be placed on the NAE, and that by allowing the surface to vary freely we are permitting unphysical solutions: for instance, allowing more material to be removed at a given time and temperature than could have cooled to that temperature given the emplacement history up to that point. Placing extra constraints on the NAE surface, or fitting a parameterized NAE surface based on physical models of lava-flow dynamics or terrestrial observations, could yield better results.

Another way of ‘regularizing’ the solution would be to integrate the prediction of the NAE, and therefore radiances, with lava-flow models themselves. The radiances are a function of the NAE, which is intimately related to the time-varying manner with which the lava is emplaced. Rather than inverting radiance for effusion rate and using effusion rate to drive lava-flow models, an iterative refinement where the flow is modelled, radiances predicted, compared to those observed, and the model refined accordingly might give a more coherent picture of the dynamics of the situation. Another issue that requires further investigation is that of temporal resolution. The degree of overlap of the kernels in the time direction increases with wavelength, such that shorter wavelengths should better resolve very rapid events such as Vulcanian eruptions.

One implication is that multispectral time series from orbital sensors can be utilized to gain an insight into the nature of the lava emplacement process, in terms of variation in rates of emplacement and removal. This depends on the availability of bands in the SWIR, MIR, TIR (and NIR if possible) with low-enough noise, sufficient spatial resolution and an appropriate dynamic range to recover the volcanic radiance time series. This presents a problem for geostationary satellites such as SEVIRI for which large footprints can result in weak volcanic TIR signals, often contaminated by strong signals from cold cloud, and NIR signals (if such a channel is present) are likely to be very weak compared to the typical range of reflected solar radiances. Volcanic radiances in the MIR dwarf those typically emitted and reflected from the surface of the Earth, so these bands frequently saturate (e.g. in SEVIRI, AVHRR). Sensors in low Earth orbit frequently carry more bands in the NIR and SWIR (e.g. MODIS, ALI, Landsat), and with specialized bands for detecting hot features in the MIR that do not saturate (e.g. MODIS); however, the lower sampling interval may not be enough to faithfully capture the variation in radiance over the course of an eruption, and the irregular sampling in space and time make constructing the kernel matrix more difficult. The use of hyperspectral SWIR and TIR sensors might substantially increase NAE resolution; however, building the kernel matrix and performing the inversion would become very computationally expensive. A framework that was capable of integrating spatially and temporally heterogeneous observations from numerous satellites in a range of bands could reveal a substantial amount about the NAE during an eruption, which could, in turn, be used to make inferences about the nature of the eruption.

However, the most immediately useful conclusions of this study regard the interpretation of TIR and TIR+MIR radiance time series of volcanic eruptions. Owing to the comparatively flat nature of the TIR kernel over the range of temperatures at which one might expect material to be emplaced and removed, the ‘smearing’ associated with the inversion conserves total NAE, resulting in a surprisingly well-recovered TAE curve – the rate of expansion in area of the hot material. This quantity, multiplied by an estimate of the material thickness, gives an estimate of the effusion rate that is robust with respect to variation in emplacement style and accounts for the cooling of previously emplaced material. A simpler approach to applying this technique in the TIR alone case might be to assume an NAE equally distributed over a reasonable range of temperatures, and calculate the average radiant cooling curve this would produce, and then simply deconvolve the radiant time series directly. Alternatively, blind deconvolution could be used to invert for the total emplacement curve and the appropriate cooling curve simultaneously. This could be applied to time series of radiance from orbital sensors, but also from terrestrial radiometers. Despite having been superseded to a degree by thermal cameras, radiometers remain a cheaper alternative and this technique could be utilized to interpret data from networks of radiometers permanently installed around volcanic fissures and flows in environments inappropriate for long-duration installation of expensive thermal cameras. It could also be used for the treatment of existing data sets.

The flatness of the kernel is due to the averaging effect of the convolution with the interpolation kernel in equation (14), which smoothes out the initial peak in radiance. The flatness of the kernel is thus a function of the cooling curve and timescale over which area emplacement varies. This implies that for a volcanic process with a given characteristic cooling curve and timescale of fluctuation, an appropriate waveband can be selected to yield a flat kernel and a robust estimate of total area emplacement. In this case, for a columnar cooling curve and a timescale of 900 s, the TIR waveband is appropriate, while the addition of another band (in this case MIR) might improve the recovery or make recovery feasible where an ideal waveband is not available for technological or atmospheric transmission reasons. The flatness across the domain is a function of the progressive overlap of kernels in the temperature direction: ironically, it is the very property that makes inversion for the NAE so difficult that facilitates the robust recovery of the TAE curve in some situations.

This research was undertaken as part of the NERC consortium project ‘How does the Earth’s crust grow at divergent plate boundaries? A unique opportunity in Afar, Ethiopia’ (grant number NE/E005535/1). CO is additionally supported by the National Centre for Earth Observation (COMET+). This paper benefited from numerous discussions with Nial Peters, and we thank Laszlo Keszthelyi and an anonymous reviewer for their helpful reviews