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

Abstract We present a novel method for interpreting time series of multispectral observations of volcanic eruptions. We show how existing models relating to radiance and area emplacement can be generalized into an integration-convolution of a Net Area Emplacement (NAE) function and a cooling function, assuming all surfaces follow the same cooling curve. The NAE describes the variation in the rate of emplacement of hot material with time and temperature, while the cooling function describes the cooling of a hot surface with time. Discretizing the integration-convolution equation yields an underdetermined matrix equation that we solve using second-order Tikhonov regularization to stabilize the solution. We test the inversion by modelling plausible NAE surfaces, calculating the radiances, adding noise and inverting for the original surface. Three or more spectral bands are required to capture the overall shape of the NAE, and recovering specific quantities is difficult. Single wavebands that yield flat kernels recover the total area emplacement curve (rate of increase of hot area – the integral of the NAE with respect to temperature) surprisingly well due to their property of conserving NAE, suggesting novel methods for calculating area emplacement rates (and effusion rates) from time series of satellite images and radiometer measurements.

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(Harris et al. , 2003Ripepe 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 twocomponent 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 lavaflow 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 finitedifference 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 highresolution 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 mm, 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.

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

Derivation
We model the radiance emitted by the lava-flow surface as the sum of the radiance of each subpixel surface element i of fractional area f i and temperature T i , radiating by the Plank function L (including emissivity, which is assumed constant for all area elements): (1) If we assume that all surface area elements follow the same cooling curve, then at time t the temperature of each surface element will be a function of the time at which it was emplaced, t i , and the radiance can be considered a function of time as well as of wavelength: The fractional area f i occupied by area element i is equal to rate of area emplacement D f i /Dt at time t i multiplied by the interval Dt i over which it was emplaced: Allowing Dt i to become very small, replacing t i with t and assuming that the eruption starts at t ¼ 0, gives: This is the convolution model, where the radiance at a time and wavelength is the convolution of the Planck function of the cooling curve with the rate of area emplacement curve. The convolution model can be considered to simulate a 'conveyor belt'-type scenario, where surface area elements are emplaced at a vent at a given temperature and then cool as they move away without changing size or being removed, as shown in Figure 1a. However, in real eruptions, not all material is emplaced at the same temperature: surfaces may be emplaced at the lava-eruption temperature near the vent, while surfaces generated at tumuli or breakouts downflow may form from cooler lava, and surfaces exposed in cracks between crustal plates will exhibit a range of temperatures from that of the cooler viscoelastic subcrustal layer to the hotter lava core. These small temperature differences are important for the thermal balance of the lava flow due to the fourth power of temperature, T, in the Stefan-Boltzmann equation. In addition, crust is removed from the field of view of a pixel by rollover at the flow front, foundering and subduction of crustal plates in lava channels, and by being overtopped by newer lava. Furthermore, the instrument observing the eruption may have an obstructed view of vents and cracks, and so fresh material may only come into view having already cooled somewhat. We can attempt to account for all these different processes by allowing surface area elements to be emplaced or removed at a number of different temperatures, T j , such that the radiance is the sum of the convolutions of the emplacement curve and cooling function for all j emplacement temperatures: where T j (t, l) and (d f /dt) j are, respectively, the temperature decay curve and the rate of area emplacement for material emplaced at temperature T j . If we assume that all surface area elements exhibit the same cooling curve irrespective of emplacement temperature, then the temperatures T j can be described as offsets in time along the same cooling curve, given by the temperature equivalent age t j .
The temperature decay function, T(t, t j ), becomes a function of two parameters -time, t, and the time offset, t j -such that T(0, t j ) ¼ T j . The area emplaced for the jth convolution curve is the rate of area emplacement integrated over both t and the temperature equivalent age interval, Dt j : As before, allowing D t j to go to 0 and replacing it with t 0 gives the equation as a double integration: The d 2 f/dt dt 0 term can be considered a Net Area Emplacement (NAE) function that gives the net emplacement or removal of the fractional surface area at time t and temperature T 0 or temperature equivalent age t 0 ; the idea is shown schematically in Figure 1b. (Note that here, as in the rest of the paper, we plot the variation in NAE against temperature rather than against temperature equivalent age (t 0 ) as, owing to the exponential relationship between T and t 0 , plotting against T compresses the function, thus making visual interpretation easier and thinking of NAE in terms of temperature is more intuitive.) Integrating over t 0 gives the rate of change of area of the whole flow at time t, which we will refer to as the 'Total Area Emplacement', or TAE, and integrating over t 0 and t gives the final total area of the flow. Multiplying the TAE by the average flow thickness (in the case of an advancing lava flow) gives an estimate of the volumetric effusion rate. Considered as a function of temperature rather than temperature-equivalent age, the NAE gives the rate of change of the temperature distribution of the lava-flow surface that is not accounted for by cooling. The NAE function for the eruption of a lava surface thus contains important information about the total area of lava emplaced, the temporal variation in the rate of lava emplacement, and the balance between emplacement and removal of the surface material at different temperatures, which in turn tells us about, for example, the role of cracking, subduction and entrainment of cold crust cooling the core in the emplacement process. The L(T(t, t 0 ), l) term is the kernel function that gives the radiance of a surface area element at time t of temperature equivalent age t 0 at wavelength l. Our assumptions about how the temperature of a surface area element cools are embodied in T(t). The challenge is to recover the unknown NAE from a time series of multispectral radiance measurements based on an assumed kernel function.
Discretizing the problem 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 temperatureequivalent 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.
We discretize wavelength by simply taking the midpoint wavelength of each band, such that for a given wavelength l: where R l (t) and L l (t, t 0 ) are the observed radiance and modelled radiance cooling curves, respectively, for a band centred at wavelength l (the temperature cooling function, T, is subsumed into L for convenience). We discretize in the t 0 direction by integrating over t 0 using simple collocation, where we sum over a number of NAE values located in a series of temperature-equivalent times t 0,i or temperatures T i : Finally, we discretize over t by replacing NAE(t, t 0 , i ), with the convolution of a linear interpolation function (a triangle function) and a Dirac comb: (10) where L(t) is the triangle function: and comb (t − t, t 0, i ) is a Dirac comb with values of the NAE to be interpolated between (the a i,j ), spaced every Dt in the t direction and at positions t 0,i in the t 0 direction: Such that the triangle function linearly interpolates between successive a i, j in the t direction. Substituting into equation (9) and rewriting in terms of the convolution operator ⊗ for simplicity gives: By the associativity and commutativity rules for the convolution operator, we set: to give: and the convolution of a Dirac comb with the kern function is just the sum of the offset kern functions weighted by the a i,j : Thus, we have simplified the function such that the radiance at a point in time is simply the sum of offset kernel functions weighted by the NAE values. If we sample the NAE surface at the same times as the radiances R l (t) are measured, we can rewrite equation (16) as the matrix equation: where R l is a vector of length M of the radiance measurements, A is vector of length N × M of the NAE coefficients a i, j arranged as a sequence of time series for increasing values of t 0 and K l is a matrix of dimensions N × M by M, containing values of the kernel function given by: The matrix K l is thus a row of N M×M Toeplitz convolution submatrices (all diagonals are the same) and only M values of k have to be calculated using equation (18) for each one, reducing the number of calculations by a factor of M. The equations for all l bands of interest can be concatenated to form a single equation: This gives the final forward equation expressing a vector of time series of radiances in all observed bands as a simple linear function of a vector describing the NAE surface: Given an assumed kernel matrix K, we aim to invert the radiance vector R for the NAE vector A.

Choosing a cooling curve
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.

Analysing the problem
We have discretized equation (7) in terms of wavelength bands and samples every Dt in the time direction and at points t 0,i in the temperatureequivalent-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, Fig. 3. Kernel vectors for a range of band wavelengths. Each plot shows the magnitude of the kernel over the NAE surface for the 50th observation in a time series of 100 observations at 15 min intervals. The left-hand column shows kernels for lava surface area elements modelled as the top of a 2 m lava column, while the right-hand column shows those for 2 cm-radius spherical clasts (parameters for both column and spherical models: eruption temperature of 1373K, ambient temperature of 300K, thermal diffusivity of 7 × 10 27 m 2 s 21 , thermal conductivity of 1.5 W m 21 K 21 , emissivity of 0.95 and convection coefficient of 60). 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 0 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 nontrivial '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.

Solving the equation
We can attempt to compensate for the ill conditioning of the kernel matrix K using regularization. This allows us to select -from amongst the broad category of NAE vectors that provide 'reasonable fits' due to the high noise sensitivity -an 'optimal' NAE vector according to some criterion. This is achieved by minimizing the sum of the square of the residuals between the predicted and observed radiance vectors, and some measure of the magnitude of the NAE vector or its slope. The idea is to minimize some measure of 'complexity' of the NAE vector such that it is just sufficiently complex to fit the data to a reasonable misfit while avoiding overfitting. We do this by solving the regularized least-squares problem given in equation (21) (Aster et al. 2011): In this study, we minimize the sum of the squares of the Laplacian of the NAE surface, which is effectively finding the NAE surface with the least total squared curvature. This is consistent with the a priori assumption that the processes determining the shape of the NAE do not change rapidly in either the time or temperature direction. Solving equation (21) subject to Laplacian regularization (also known as second-order Tikhonov regularization: Aster et al. 2011) involves solving the following least-squares equation: where L is the discrete Laplace operator for the 2D NAE surface defined by the vector A, C −1 N is the inverse of the radiance covariance matrix (a diagonal matrix with the noise variance for each band on the leading diagonal and zeros off the diagonal), and a is the weighting parameter that determines the relative contribution of the data misfit and model Laplacian semi-norm to the cost function. We solve this equation using the QR decomposition approach. The value of a is chosen according to the L-curve criterion, where equation (22) is solved for a range of values of a, each of which specifies an optimal ||K A − R|| 2 2 and ||L A|| 2 2 . Plotting the log of these two values against each other usually gives an L-shaped curve, where the corner of the L specifies the value of a above which further 'simplification' of the solution model leads to a rapid increase in misfit (Aster et al. 2011). We select this point by solving equation (22) for a range of values of a, interpolating between the corresponding (log||K A − R|| 2 2 , log ||L A|| 2 2 ) pairs, finding the point of maximum curvature, finding the value of a at that point and solving equation (23) with that value. An example L curve is shown in Figure 4.

Test models
We can test the degree to which we can recover the NAE surface using the aforementioned regularization scheme by constructing plausible NAE surfaces, calculating the expected radiance, adding simulated sensor noise, and then solving equation (22) and comparing the original and recovered models. In the absence of observed NAE surfaces or parameterized surfaces derived from numerical models of lava emplacement, we define a plausible functional form a priori. We assume that the NAE surface varies smoothly in time and temperatureequivalent age since there are a wide range of different processes acting to emplace and obscure surfaces at different temperatures and different geomorphologies during an eruption (see Figure 1), yet we also assume that the NAE is localized since the NAE before and after the eruption, and for temperatures above the eruption temperature and below the ambient temperature, must approach 0 rapidly. These requirements are satisfied by the product of two skew-Gaussian functions: yielding a smoothly varying but localized 2D surface defined by six parameters specifying its position on the NAE surface (m t , m t 0 ) width in both directions (s t , s t 0 ) and skew in both directions (b t , b t 0 ) The skew parameter is particularly useful as it can be used to approximate the asymmetric waxing and waning pattern often seen in effusion rates from fissural eruptions (Wadge 1981). The sum of a number of these pairs can be used to build up an NAE surface consisting of positive regions of net area emplacement and negative regions of net area removal:

Results
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 0 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 27 m 2 s 21 , thermal conductivity of 1.5 W m 21 K 21 , 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 21 s 21 . Radiances are then calculated from these NAE surfaces using equation (20), and Gaussian noise  (24). The radiance time series is then calculated using equation (20). (b) L-curve found from least-squares fits to the radiances by solving equation (22) for a range of regularization parameters. The regularization parameter associated with the point of maximum curvature of the L-curve is found (red dot) and the least-squares solution is calculated for that value, giving the regularized 'best-fit' solution. (c) The radiances calculated from the test NAE model plus added noise (blue curves) and the fitted radiances from the retrieved NAE surface (red curves). (d) The retrieved NAE surface. 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 mm, 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 mm) 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 0 /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.

Discussion
The test models presented do not exhaustively sample the full space of possible NAE surfaces; The total area emplacement, found by integrating the NAE over temperature. Lower panels: recovered NAE surfaces and total area emplacement curves for inversions using various combinations of bands (S, shortwave; M, mid-infrared; T, thermal infrared; N, near-infrared). The solid curve shows the model total area emplacement; the dashed curve shows the recovered total area emplacement. NAE surface and total area emplacement plots have the same axes as the plot for the original model, omitted here for clarity. 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, Fig. 6. A simple scenario with a positive NAE surface at moderate temperatures and no removal of material, with the columnar cooling kernel. (a) The modelled NAE surface. (b) The total area emplacement, found by integrating the NAE over temperature. Lower panels: recovered NAE surfaces and total area emplacement curves for inversions using various combinations of bands (S, shortwave; M, mid-infrared; T, thermal infrared; N, near-infrared). The solid curve shows the model total area emplacement; the dashed curve shows the recovered total area emplacement. NAE surface and total area emplacement plots have the same axes as the plot for the original model, omitted here for clarity. 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%

Conclusions and implications
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) The total area emplacement, found by integrating the NAE over temperature. Lower panels: recovered NAE surfaces and total area emplacement curves for inversions using various combinations of bands (S, shortwave; M, mid-infrared; T, thermal infrared; N, near-infrared). The solid curve shows the model total area emplacement; the dashed curve shows the recovered total area emplacement. NAE surface and total area emplacement plots have the same axes as the plot for the original model, omitted here for clarity.
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 The total area emplacement, found by integrating the NAE over temperature. Lower panels: recovered NAE surfaces and total area emplacement curves for inversions using various combinations of bands (S, shortwave; M, mid-infrared; T, thermal infrared; N, near-infrared). The solid curve shows the model total area emplacement; the dashed curve shows the recovered total area emplacement. NAE surface and total area emplacement plots have the same axes as the plot for the original model, omitted here for clarity. 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. Fig. 11. A complex scenario with a positive NAE surface at high temperatures and substantial removal of material at lower temperatures, with the columnar cooling kernel. (a) The modelled NAE surface. (b) The total area emplacement, found by integrating the NAE over temperature. Lower panels: recovered NAE surfaces and total area emplacement curves for inversions using various combinations of bands (S, shortwave; M, mid-infrared; T, thermal infrared; N, near-infrared). The solid curve shows the model total area emplacement; the dashed curve shows the recovered total area emplacement. NAE surface and total area emplacement plots have the same axes as the plot for the original model, omitted here for clarity.
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.