Atmospheric Measurement Techniques Volcanic SO 2 and SiF 4 visualization using 2-D thermal emission spectroscopy – Part 2 : Wind propagation and emission rates

A technique for measuring two-dimensional (2-D) plumes of volcanic gases with thermal emission spectroscopy was described in Part 1 by Stremme et al. (2012a). In that paper the instrumental aspects as well as retrieval strategies for obtaining the slant column images of SO 2 and SiF4, as well as animations of particular events observed at the Popocat́ epetl volcano, were presented. This work focuses on the procedures for determining the propagation speed of the gases and estimating an emission rate from the given image sequences. A 2-D column density distribution of a volcanic gas, available as time-consecutive frames, provides information of a projected wind field and the average velocity at which the volcanic plume is propagating. This information is valuable since the largest uncertainties when calculating emission rates of the gases using remote sensing techniques arise from propagation velocities which are often inadequately assumed. The presented reconstruction method solves the equation of continuity as an ill-posed problem using mainly a Tikhonov-like regularisation. It is observed from the available data sets that if the main direction of propagation is perpendicular to the line-of-sight, the algorithm works well for SO2, which has the strongest signals, and also for SiF4 in some favourable cases. Due to the similarity of the algorithm used here with the reconstruction methods used for profile retrievals based on optimal estimation theory, diagnostic tools like the averaging kernels can be calculated in an analogous manner and the information can be quantified as degrees of freedom. Thus, it is shown that the combination of wind field and column distribution of the gas plume can provide the emission rate of the volcano both during day and night.


Introduction
Measuring volcanic gas emissions is important for several reasons.Gas fluxes are commonly monitored by governmental agencies and research institutions for disaster prevention purposes, especially in the case of active volcanoes near cities and along major flight routes.Also, a detailed and current knowledge of the emissions of volcanoes, both eruptive and passive degassing, is important as it is well known that SO 2 , for example, has an indirect impact on radiative forcing and is therefore needed in climate model studies on regional and global scales (Halmer, 2002;Robock, 2000).
However, only a small number of active volcanoes are studied worldwide and the frequency of measurements is strongly dependent on the method used and the capacities and funding available by the responsible local institutions.Nowadays, much can be done using satellite remote sensing, but these measurements have usually large uncertainties and need continuous improvements in the quantification of gas emissions, preferably through validations with ground studies.The most monitored volcanic gas is SO 2 due to its abundance and the fact that it can be readily measured using common analytical methods like DOAS (Differential Optical Absorption Spectroscopy) in the ultraviolet or FTIR (Fourier Transform Infrared Spectroscopy) in the infrared, both in absorption or thermal emission modes.
Each technique and measurement geometry needs a different treatment to derive information about the emissions from the measured quantities.For example, a nadir observation from space can obtain a snapshot of the horizontal distribution of SO 2 columns, and with the use of an inverse model the emission can be estimated (Kristiansen et al., 2010).Near the ground, a traverse around the volcano using scattered sunlight and the DOAS technique from various platforms gives a cross-section of the volcanic plume as a result.This can be done, for example, by walking (McGonigle et al., 2002), from an ultra-light aircraft (Grutter et al., 2008a) or by using a fixed instrument to scan the plume from below, as was done by the NOVAC project in a network of multi-axis instruments (MAX-DOAS) set up in various volcanoes around the world (Galle et al., 2010).
However, apart from the gas column measurements, which represent the plume cross-section, the wind velocity also has to be known to calculate emission rates.Normally this information is taken from ground measurements, radiosondes or reanalized meteorological data products from models.Correlations between measurements from three instruments on different locations (McGonigle, 2005), or measuring at different angles (McGonigle et al., 2005) has been done to determine the velocity.Recently, Johansson et al. (2009a) developed a dual-beam device which allows the estimation of the wind velocity with a single device.
Instrumentation for imaging gas distributions in volcanic plumes is a current topic, since it aims to provide information about the dynamics and fluxes as has been done by means of UV cameras (Bluth et al., 2007;Kantzas et al., 2010).Most scanning and imaging devices use UV radiation; however, spectrally resolved infrared radiation also has been used for characterizing the spatial distribution of industrial and volcanic emissions (Harig et al., 2005;Grutter et al., 2008b;Gross et al., 2010;Stremme et al., 2012a).
In this contribution, it is shown how the propagation velocity of the gas plume can be determined from a sequence of two images and used to calculate an emission rate around the volcano.Column measurements have been retrieved from thermal infrared spectra taken with a FTIR spectrometer at 4 cm −1 resolution.A scanning device has allowed for the recording images of the slant column density of SO 2 and in some favourable cases also of SiF 4 (see Part 1, Stremme et al., 2012a).Figure 1 shows an example of two sequential images with SO 2 columns in false colours from which the emission rates have been calculated.
This work presents a new method which reconstructs the 2-D projection of the trace gas flow and its sources simultaneously by the exploitation of the continuity equation.The method uses a simple physical forward model based on the continuity equation and an inversion based on Thikonovregularisation and optimal estimation (Rodgers, 2000).The reconstruction of the plume propagation and sources does principally not depend on any external input about wind direction, but its formalism allows for the use of external information in addition, as shown in this work.
One advantage of the use of the continuity equation is that a source and sink distribution is being retrieved simultaneously.The amount of information gained by this scheme can be chosen by adjusting the strength of the regularisation, and the retrieval strategy can be adjusted to the specific application by taking into account frame rate, signal/noise ratio of the measurement and the meteorological conditions.The analysis of the recorded volcanic gas images uses a multi-step retrieval strategy.The strategy developed here uses an additional correlation pattern to optimise the retrieval for the calculation of the emission rate.

Measurements and instrumentation
The measurements considered here were taken with a Scanning Infrared Gas Imaging System (SIGIS).It is composed of a two-dimensional scanning mirror coupled to a video camera.The mirror sends the thermal radiation to a telescope and an FTIR spectrometer (OPAG 22,Bruker Daltonics,Leipzig,Germany).The instrument uses an MCT detector with a stirling-cooler and covers a spectral range from 600-5000 cm −1 .Its best spectral resolution of 0.5 cm −1 has been used for solar absorption measurements and for thermal emission spectroscopy (Grutter et al., 2008a,b;Stremme et al., 2009Stremme et al., , 2011Stremme et al., , 2012b)).However, for imaging and the analysis considered here, the spectral resolution of 4 cm −1 is more convenient.The measurements were performed from the Altzomoni Atmospheric Observatory, which is around 12 km north of Popocatépetl at an altitude of about 4000 m a.s.l.Slant columns of SO 2 and SiF 4 referring to an almost horizontal line-of-sight were retrieved from the IR spectra.For more details of the instrumentation, used algorithm and realisation of measurements, see Part 1 of this work (Stremme et al., 2012a).

Methodology
The spatial distribution provided by the column density images together with their temporal change allow for estimating the propagation speed and direction as well as the source  strength of the plume.The approach for how this information can be extracted mathematically is the main topic of this work.

Overview
A simple, linear and physical forward model based on the continuity equation is used to explain two consecutive images.The use of a physical forward model is the main difference to other methods that use correlations (e.g.Johansson et al., 2009a) and/or external information as, for example, the wind information from a weather forecast model (e.g.Lübcke et al., 2012).The result allows estimation of the trace gas flux (slant column × velocity) perpendicular to the lineof-sight (LoS) in each pixel.From two consecutive images a distribution of the emission rate of the volcano as a function of time can be reconstructed.
Our retrieval method seeks the reconstruction of an optimised field of 2-D wind components and a field of sources, but the information available from the measurements is limited.According to inversion theory, the retrieval can be optimised for different purposes (e.g.Steck and v. Clarmann, 2001) and the sought parameter is not always optimally retrieved.We use a multi-step retrieval, whose steps differ with respect to the input and constraints used but all rely on the same basic principle, which is the inversion of the continuity equation described in Sects.3.2, 3.3, and 3.4.In the first retrieval step, a 2-D field of wind velocities and source strengths is determined from the measured slant path column amounts.From the calculated wind velocities, the mean wind direction is inferred, which is used to constrain the second retrieval of wind velocities and source strengths.From these, emission rates at the crater are inferred, and by application of a cross correlation method (Sect.3.5) the mean windspeed is inferred, which is used to constrain the final retrieval of wind velocities and source strengths from which the final emission rates are inferred.This sequence of operations is summarised in Table 1.More insight why a multi-step retrieval scheme is used is given in Sects.3.5 and 4, where a description of how the emission rate is obtained from the retrieval result is also included.
With "emission rate" we refer to the emission of the volcano, which results from the trace gas flux integrated over specific plume cross-sections.These can be compared with the integrated sources inside the images while dropping the sinks near the frame borders.

Basic principle
The density distribution and its temporal evolution is described by the equation of continuity.This applies not only to the number density (ρ) of SO 2 , SiF 4 or any detectable gas with volcanic origin, but also to the slant column densities (cl:= LoS ρ • ds) measured by the instrument.The ds is the infinitesimal path along the LoS of the instrument to the top of the atmosphere, to a sight-limiting cloud or to the body of the volcano.The continuity equation in its general form reads Integrating Eq. (1) along the line-of-sight yields Since there is no hope for retrieving the wind and source component along the line-of-sight, we only consider the two dimensions (x, y) perpendicular to the line-of-sight.Disregarding mixing processes, we replace the respective wind components (v x , v y ) by their line-of-sight averages weighted by molecule number density, where the ∇-operator and the velocity v 2-D are reduced to two dimensions,  (b) is the same as (a) but calculated using the forward model K and the solution vector x, consisting of the retrieved vector field and the retrieved source field.The residual (c) is the difference between the plots above (observed -simulated).Due to the choice of S −1 e (see text, Sect.3.4), the residual of the column differences are larger outside of the plume.The location of the Popocatépetl volcano is indicated.more details can be found in Appendix A.
For practical purposes, the 2-D approximation has a minor impact on our retrievals as long as number concentrations of SO 2 (or SiF 4 ) outside the plume are negligible.
A two-dimensional wind field (v 2-D ) might be reconstructed by solving the continuity equation for the column densities of a gas with volcanic origin using Eq.(3).
The left term of the equation represents the change from one to the next 2-D distribution of the observed gas column density (cl) per time interval, as is shown in Fig. 2a.v 2-D is the 2-D propagation vector and q is the source term that describes the appearance of new gas signals in each pixel of the image, being equivalent to LoS ∂ ∂t ρ • ds in Eq. ( 2).The source term has the unit of molecules per area per time, which is equivalent to emission flux.The term accounts also for chemical processes as, e.g.SO 2 destruction through photolysis or dissolution into the aqueous phase.

Forward model
To apply the physical Eq.(3) to measured data, discretisation in time and space are necessary.The smallest available time interval dt is given by the frame rate (1/t frame ) and the smallest distance in the x-y direction (where x represents the Fig. 3. Plume propagation vectors obtained from two sequential gas column images.SO 2 columns are retrieved from measured IR spectra assuming a temperature of 269 K at 5500 m a.s.l.obtained from a radiosonde at 06:00 LST.The wind field is retrieved from this and the following frame, solving the equation of continuity.The dark solid line represents a linear approximation for the trajectory and the cross-sections perpendicular to it.The cross-sections are separated by one-minute intervals and the velocity of the propagating plume is projected to the drawn trajectory. horizontal and y the vertical direction) is given by the angle that the motor moved the mirror from one to the next pixel of the 2-D image, multiplied by the distance to the plume.The pixels are quadratic with a size x = y of around 120 m at a estimated distance of 12 km, and their positions are described by the subscript (i,j ) where the position in the raw is given by the first index i and the position in the column by the second index j .dcl dt if both neighbours are available if both neighbours are available. (5) At pixels at the edge, asymmetric finite difference quotients are used instead.Equation ( 6) represents a linear forward model where y is the measured quantity, x is the solution vector and K the linear forward model.
For images with n × m pixels, Eq. (3) describes after its discretisation in Eq. ( 6) an under-determined problem with only n × m equations, corresponding to the differences in the columns of two sequential frames y = (dcl 11 , . . ., dcl mn ) and the three times larger solution vector x = (v x 11 , . . ., v x mn , v y 11 , . . ., v y mn , q 1 , . . ., q mn ) containing the wind components (v x ij , v y ij ) and the m × n source strengths (q ij ) (cf.Figs. 3 and 4b,respectively).
The left and right side of Eq. ( 6), shown as Fig. 2a, b, represent the measured and simulated differences of the As the units are arbitrary, they have to be compared with the constraint (R) concerning the wind components and with the weighting of the measurement (S −1 e ).(b) Retrieved sources are also part of the solution vector.As the wind field is constrained, this additional fit parameter improves the result significantly.The clean atmosphere is constrained with a smoothing constraint, while the pixels of the borders and the volcanic body are constrained with an optimal-estimation-like regularisation (R-matrix is diagonal).In this example, 11.7 independent pieces of information are included in the vector, which has 714 elements to describe the sources.As expected, the retrieval finds the strongest sources above the crater.column densities between two consecutive images, respectively.As the propagation during the time between two images (dt) might be larger than the pixel size, the images can be smoothed before dcl dt and ∇cl are calculated.With this, not only neighbouring pixels are taken into account but the smoothed gradient (∇cl) might be more adequate for the finite time interval between two measurements.
The matrix It is rather sparse and given by the following expression: The location of the element (v x 2-D ) l is found as (v x ij ) in the image using the 2 indices (i,j ) by the relation l(i, j 1 Eqs.( 8) and ( 9) are not valid for pixels (i,j ) at the edges of the image, as the form of the gradient is different.
The first expressions in K x and K y are diagonal matrices, but the second term expresses the partial derivation of x and y the will operate on the 2-dimensional wind field.

Inversion
As the solution vector has a larger dimension than the measurement vector, the problem would be ill-posed even if there were no measurement errors.The ill-posed Eq. ( 6) is solved by constrained least squares fitting with the following penalty function (Eq.11) to be minimized: x = x ret is the solution which minimizes the penalty function in Eq. ( 11) and contains the retrieved wind propagation vectors and the source strength in all pixels.
The covariance matrix S e −1 characterises the uncertainties of the slant path column amounts and is important for the optimal weighting of measurements in Eq. ( 11); for those pixels where the column density errors are larger, the residual element (y − Kx) i,j has less impact on PF (cf.Fig. 2).In this work, S e −1 is chosen to be diagonal, which means that the errors of slant path column measurements are uncorrelated.The diagonal elements of weighting matrix S e −1 are estimated by the following ad hoc approach: The Pearson's correlation coefficient is calculated between the measured spectra from which all fitted interference gases, offset, slope and curvature have been subtracted and a reference spectrum (see Part 1, Stremme et al., 2012a).The correlation coefficient as a data matrix represented in a false-colour image is a possible output of the used algorithm and therefore the basis for the weighting matrix.The coefficient r 2 Pearson -together with a user-defined threshold (r 2 threshold ), r 2 max the highest correlation coefficient in the frame and a simple function ( through a smoothing constraint (Tihkonov-type regularisation), (ii) adding a priori information (Bayesian approach) and (iii) parametrisation of the solution vector (Tikhonov, 1963;Rodgers, 2000).The construction of the matrix R takes aspects of these three into account and the borders between the strategies are rather blurred in an empirically adjusted retrieval.
The first order difference operator L 1 (e.g.Steck, 2002;von Clarmann and Grabowski, 2007) calculates the differences between neighbouring values.λ • L T 1 L 1 constrains the solution towards a smoothed solution without constraining the average values of the x-,y-components of the velocities and sources.The optimal estimation strategy is realised partly and in an empirical manner by another term which constrains (with strength µ) each value of the solution vector towards an a priori value of x a .This constraint affects the average wind velocity of all pixels, but it allows inclusion in a simple manner a priori information of independent measurements or even from the same data by applying a different strategy like a cross-correlation (see Sect. 3.5).In the algorithm used in this work, only one a priori velocity (v x a , v y a ) valid for all pixels in the image is chosen.
Information concerning the location of sources and sinks is also needed.In the solution of Eq. ( 3) no limits of the area are taken into account; however, at the border of the considered area a SO 2 plume is expected to enter or leave the observed area.This should be reflected by the retrieved sources or sinks in exactly those pixels where the plume meets the border.The constraint should ensure that such a solution is not suppressed.Therefore, a term is added to free the border and all pixels which are considered as potential sources from a smoothness constraint and treat their regulation individually (see Fig. 4a).The choice of location where potential sources are expected and constraining differently as background might be considered a kind of parametrisation.R is constructed from three terms each consisting of three block matrices: Two blocks concern the x and y components of the wind field (v x 2-D , v y 2-D ) and the third the sources q.
The first term of R is a smoothing constraint which encourages the neighbouring wind velocity components and sources in x to have a similar size and therefore the wind vectors to point in the same direction.The strength of the constraint for the wind field is given by λ x , λ y , which are chosen to have the same strength.The strength of the field of sources λ src is independently chosen but has the same smoothing effect.
The second term constrains the solution towards the a priori.This constraint is first (as here in the first two retrievals) turned off (λ x,y = 8000, µ x,y = 0) for the wind velocity components, but it offers the possibility to include other information about the wind velocities taken from a crosscorrelation, as in the example shown in this work, from the average of prior retrieved velocities or from independent measurements, e.g. from radiosondes.Here the final retrieval is realised without smoothness constraint (λ = 0) and constrained towards the wind velocity which has been calculated from a cross-correlation (λ x,y = 0, µ = 100).The sources are constrained slightly towards zero.
The third term of the constraint refers to the sources only.The matrix R src is automatically initialized using the information where sources or sinks are expected.The whole frame is a potential source or sink as the SO 2 plume can enter or leave the image.The same argument applies to the body and crater of the volcano and opaque clouds.Both can be identified from the recorded spectra as the observed intensity in the mid-infrared is much higher in these cases.For these pixels the constraint in terms 1 and 2 will be eliminated and just a very small value is added to the diagonal to ensure that the problem can be solved mathematically.

Calculation of the emission rate and cross-correlation
Traditional methods use a cross-correlation approach to obtain an average propagation vector (Johansson et al., 2009a;McGonigle et al., 2005).In our approach, we retrieve the individual 2-D velocity vectors for each pixel but the crosscorrelation method still is useful for correction of the average wind speed.Usually two time series measured at different locations are used for the cross-correlation.In this work two time series are reconstructed with the help of a 2-D wind field from spatial distribution recorded at two times.First, the 2-D wind field is retrieved using a pure smoothing constraint and the spatially averaged wind velocities are calculated (for the average the wind is weighted by the SO 2 column).Typically the direction is well retrieved but the absolute wind speed might be off.After the direction is retrieved, the upper boundary of wind speed is constrained to a moderate value of 2 m s −1 and a new wind field is obtained (arrows like in Fig. 3)2 .
The SO 2 flow column density (j ) is calculated by multiplying the wind field with the SO 2 column densities in each pixel.As a sequence of two images are used for the 2-D wind field retrieval, either the first (former) or the second (latter) image and 2-D column distribution cl former , cl latter can be used for an emission rate calculation.) at different distances from the crater represented by the cross-sections drawn in Fig. 3 and using the reconstructed wind field.The distance is shown as the time of its propagation after the starting point of the linear trajectory (Fig. 3).The solid trace shows the flux calulated using the first image of SO 2 columns and the dashed line is the flux at the same time, but based on the next SO 2 image.The reconstruction of the wind field is done using a constraint which limits its spatial resolution.
Especially near the crater, the dispersion of the SO 2 plume is not very well represented by the retrieved wind field and the SO 2 emission is underestimated (more or less below 6 min).A similar effect occurs at the other end of the image.In the centre of the plume, where the linear trajectory is well represented, both SO 2 emission estimations are consistent and oscillate around 1000 t day −1 in this case.
Below is the same equation but in other coordinates; k is along the trajectory and s is perpendicular to it.
According to the law of Gauss, the source strength or emission rate is the total flux ρv through a closed surface.The virtual surface (e.g. a cone with the observer in the top, or a cylinder or prism for large distances between observer and plume) should enclose the source, but only the part in which the surface intersects the plume contributes to the integral.This law simplifies further when working with column densities, which means that one integration over the 2-D surface is already realised and only the integration along the baseline located in the image plane is missing.On the basis of this principle, vertical column measurements on closed or downwind traverses around industrial complexes (e.g.Rivera et al., 2009), mega cities (e.g.Johansson et al., 2009b) and volcanoes (e.g.Grutter et al., 2008a) have allowed for the determination of the total emission rates.The same principle applies if a scanning device measures a closed traverse of almost horizontal slant columns around a point source.An image of slant columns allows choice of different cross-sections and each traverse can be parametrised due to the estimated time-distance from the carter (see Eq. 19).
A time series of the emission at different distances from the crater results from integrating the flux densities perpendicular to the trajectory (near-horizontal line in Fig. 3).This calculation can be done separately for the two false-colour images (former and latter), so that two time series of emission can be calculated.
n is a normalised vector with direction of the propagation of the plume and ds k as a small part of the cross-section (C k ), which is perpendicular to n and has the distance to the crater d k .The location of the crater is the beginning of the linear track (black, straight line in Fig. 3).The corresponding time t k is obtained from the distance d k and the averaged wind speed |v ret | between the centre of the source (crater), the estimated location of the crater and the beginning of the linear trajectory, and the averaged wind speed |v ret |.
The two slightly different time series are compared in Fig. 5.One uses the column density from the first frame measured and the second uses the second, which is needed for the wind field calculation.For testing the calculated wind fields, the time series are cross-correlated (Fig. 6) and the resulting time shift (t shift ) should match the time difference between the measurements (t frames ), as shown in Fig. 6.If these times are not equal, a scaling factor for the wind speeds is calculated and used for the calculation of the a priori information (v a ) in the final (third) retrieval.

Diagnostics and errors
Independently of how the wind velocities are calculated, three problems might arise if conditions are not favourable: (a) the wind speed is too fast for the frame rate and size of the recorded imagines; (b) there is a repetitive pattern with more similar maxima or minima, so that the algorithm might not distinguish them; and/or (c) there might not be sufficient structure in the images.
For case (a) the wind-speed range is limited and thus this criterion should be taken into account when choosing the frame rate while measuring, as frames need to be refreshed soon enough to capture the wind propagation features.However, if more puffs are available in the images (case b) so that the propagation of the puffs are ambiguous, different approaches could be tried.

A. Krueger et al.: Thermal emission spectroscopy of volcanic gases -Part 2
cross/correlation Time (minutes) Fig. 6.The cross-correlation of the emission time series in Fig. 5 shows a maximum corresponding to the delay between both frames.If the wind field is well retrieved, the delay obtained from this crosscorrelation (diamond at the maximum after 2.08 min) should coincide with the known time delay between both measured frames (vertical line at 2.07 min).Normally, a small difference could always be expected since the wind field result is smoothed and it is assumed that the plume does not change its direction of propagation, but larger discrepancies in the delays would mean that the wind fields have not been properly retrieved.
First, an adequate smoothing might eliminate smaller minima and maxima, which would allow the main structure to remain.If, however, a new puff appears near the crater with a similar magnitude in its anomaly as the one already in the image, the stronger smoothing might not work.In these cases a wind velocity might be obtained which would point in the wrong direction.Since such an event might appear just once in a while where most pairs of images lead to a correct velocity, an iterative analysis of the time series or Kalmannfiltering might be used to eliminate results and improve the wind speed retrieval in the time series.
In case (c), which describes the situation where there might be no sufficient structure in the plume to retrieve a wind velocity field, there is still a chance to reconstruct a mean velocity that includes the main direction and speed of propagation.This result can be later checked or improved by a cross-correlation (Sect.3.5); how much information about the wind velocity field is extracted in a retrieval is described in Sect.4.1.

Averaging kernel and resolution
The information obtained from the analysis depends on the trace gas and if the velocity is somehow perpendicular to the LoS.Therefore, the conditions can be more or less favourable for the wind field retrieval.The averaging kernel (A), also called the data resolution matrix in other research fields, is a useful diagnostic tool to track the origin of the information in a constrained retrieval.It can also be used to quantify the information in the retrieval.
where (ij ) and (kl) are indices with values from 1 to (3 × n × m) .The A does not have a block-matrix form, but here just the 3 blocks concerning the x-and y-components of the velocities and the source strength of the pixels are considered separately.Further, the discussion is limited to the diagonal of the A so that we might address the diagonal elements of the averaging kernel of the pixel i,j as A ij x ,A ij y and A ij src in an intuitive way and later at the end of the section we will briefly discuss interference between the different retrieved quantities.
The diagonal of the A reflects the sensitivity of the retrieval (Funke et al., 2009).The sum of the diagonal, trace (A), gives the degrees of freedom (DOF) of the retrieval (Rodgers, 2000) and its value shows how many independent quantities are retrieved.DOF, calculated separately for the blocks of the x-,y-components of the velocity and the source strength, show that the information is not dispersed homogeneously over the three blocks.The DOF x , DOF y and DOF src describe the independent quantities and have the values 68.6, 92.7, and 11.7, respectively.Different values are expected for the sources as they have a different character and are differently constrained.However, a difference between the DOF regarding the v x and v y components is of outermost importance.It might explain why the retrieval works well for the calculation of the wind direction, although in some observed cases there are difficulties for the proper determination of the wind speeds.
The forward model K consists mainly of the gradient of the column densities, and especially in the case of a continuous emission, this gradient is strong perpendicular to the propagation but small or even not significant in its direction.
In the example shown in Fig. 3, the main direction of propagation is in the direction of x, explaining the imbalance in the calculated DOFs for the x v x and x v y components.
The information in the retrieval concerning A x is smaller if the gas is emitted continuously and there is only a small gradient of column density in the direction of the plume.The cross-correlation described in Sect.3.5 is useful in such cases since a reliable mean velocity can be calculated and used as a priori in the subsequent retrieval.
The averaging kernel regarding the source strength retrieval is different.In the free atmosphere this term might just help to improve the fit as no other sources are expected for SO 2 .However, where the plume leaves the image a strong sink should be retrieved.The diagonal of the averaging kernel regarding the sources is shown in Fig. 7c and reflects that the retrieval is sensitive in the pixel where the S −1 e allows taking information from the images and for which the sources are not constrained by the regularisation R. Both conditions are met where the plume leaves the image.
The used regularisation regarding the strength of the sources in each pixel is visualised as the corresponding part of the diagonal of the matrix R in Fig. 4a.This figure The first two retrievals with stronger Tikhonov and with negligible OET (diagonal) constraint, (II.)Final retrieval with negligible Thikonov constraint but stronger OET (diagonal) constraint.As the main wind direction is along the x-axis, A x represents the sensitivity to a deviation from the wind speed (|v| − |v a |).As the main wind direction is along the x-axis, A y represents the sensitivity to a deviation from the wind direction (dv y = |v a |dφ).The sensitivity to sources depends on the constraint (R) and the weighting of the measurements (S −1 e ).
illustrates how the information about the body of the volcano, taken from the same frame of spectra, is used in the inversion.

The spatial distribution of source strength
The retrieved dispersion of source in Fig. 4b is part of the retrieved vector.In Fig. 4b the source dispersion shows a maximum above the crater and a strong sink, where the plume leaves the image.The averaging kernel diagonal, Fig. 7c, shows sensitivity of the retrieval in those pixels where the sink is retrieved.The maximum of the source strength above the crater is in a region where the sources are constrained by a smoothing constraint and therefore smoothed out.The location of the source is the beginning of the plume and consistent with the implemented forward model.The constraint suppresses solutions where the direction of propagation changes strongly.Therefore, the retrieved solution does not show that the SO 2 leaves the crater, but finds that a source above the crater is releasing SO 2 .The retrieved sources along the plume decrease, which might reflect that the velocity of the SO 2 molecules along the plume increases, which means that the true velocity is overestimated near the crater and near the right border it is underestimated.This might result from the regularisation of the velocity towards a mean velocity.

Errors of the emission rate
Analogously to atmospheric profile retrievals and according to Rodgers (2000), three different kind of errors should be discussed: (a) the forward model (parameter) errors, (b) the measurement-noise error (random error), and (c) the smoothing error, which mainly describes the effect of the constraint.In addition, we briefly address (d) the error due to interference between the different retrieved quantities (e.g.Sussmann and Borsdorff, 2007).

Errors regarding the forward model
The error in the emission rate can either originate from the column retrieval of the target gas or from the determination of the velocity and size of the plume.The errors affecting the column retrievals arising from the radiometric calibration and the uncertainty in the estimation of the plume temperature.Both are already discussed in detail in Part 1 (Stremme et al., 2012a), where the resulting error in the column densities were estimated to be around 20 %.Important to consider here is a geometrical error, which arises from the uncertainty in the distance between the observation site and the plume.The error (δr = r esitmated − r true ) in the estimation of the distance r affects the estimation of the amount of a gas cloud by (δr/r) 2 (see Part 1).Equally, the movement of the gas cloud can be erroneously calculated if the estimated distance r is wrong.In fact, only the change of the angle in which the plume appears dφ dt can be obtained from measurement (Fig. 6) and has to be multiplied by an estimated distance r.
The apparent velocity results to be different from the true plume velocity, and the difference is given by dφ dt • δr.The uncertainty of the emission rate thus depends quadratically (δr/r) 2 on the uncertainty in the instrument-plume distance r.A 10 % overestimation of r (δr ≈ +1 km) results in an emission rate overestimated by 21 %.No error arises from the projection, as this projection of the wind velocity is compensated due to the use of a slant slice of the plume cross-section instead of a slice perpendicular to the main wind direction.By a wind direction of 45 • towards the observer, the slant cross-section will be 1 cos(45 • ) = √ 2 larger than the cross section perpendicular to the direction of propagation, but the wind speed projected to the 2-D plane will be 1 √ 2 = cos(45 • ) smaller.However, a change in wind direction in between the two scans is not described by the forward model and will result in an error which is increased if the wind direction is almost parallel to the line of sight.Rodgers (2000) distinguishes the errors regarding the forward model as "forward parameter errors" and "forward model errors".The forward model assumes that the images of the column densities are snapshots, but the instrument is continuously scanning, so that the measurements of the pixels in one image do not belong to exactly the same time.This is, according to Rodgers (2000), a real "forward model error" and can only be evaluated if an improved forward model would be available; the errors discussed above are classified as "forward parameter errors".

Measurement-noise error
Here both the input vector dcl and also the K-matrix, which contains the gradient of the measured cl, are affected by the measurement noise (the random errors in the retrieved columns).Nevertheless, the error can be evaluated by adding artificial noise to the images and comparing the results with and without adding noise.The 2-D noise patterns for the two images are not calculated from the same but from the next following two slant column images by subtracting the smoothed images (average of the values in each pixel and in the four neighbouring pixels) from the original images.Four images are involved in the exercise, two images are used for the retrieval and two images are used to calculate independent noise for the other two images.The root mean square obtained from the noise pattern in the relevant pixels results in 13 % of the mean relevant SO 2 values (relevant refers to pixels with correlation coefficients above the threshold, Sect.3.4).The resulting error in the final retrieval was estimated from an ensemble of retrievals, for which the two noise patterns have been added and/or subtracted from the two original images and the resulting mean emissions rate observed.A typical change of the emission rate by 20 % was found and indicates that the measurement error is one of the largest error sources.

Smoothing error
As an averaging kernel, A is calculated in Sect.4.1, and a "general" smoothing error SM could be calculated, if (i) the retrieval had been constrained towards a climatological mean in the sense of optimal estimation, if (ii) the a priori covariance of the wind field and source distribution S a would be known, and (iii) the averaging kernel would be constant.However, no S a is known and the A is variable so that the discussion is limited to qualitative aspects and how the smoothing error affects the retrieved wind direction and speed, as shown in Eq. ( 23): It is clear that the wind field in Figs.8a, b is constrained and that a true velocity field would show a change in the direction with altitude above the crater.The effect of a constrained wind field is visible in Fig. 5, as the wind speed near the crater is overestimated while near the border seems to be underestimated.Such a solution results in less gas near the crater and more gas near the borders, if this residual would not be compensated by the field of sources, as shown in Fig. 4b.The smoothing errors on the emission rate times series could be calculated using Eqs.( 23), ( 15) and ( 17), if the x true could be estimated.

Interference error and dependence
Forward model parameter error (Rodgers, 2000) concerning the interference quantity and too strong constraint of the interference quantity (Sussmann and Borsdorff, 2007) can lead to errors in the target quantity of a retrieval.While the first type of interference error is commonly associated with inaccurate spectral line parameter or spectral line models (e.g.how water vapour affects methane column retrievals), the later mentioned interference error type relevant is if the retrieval reconstructs distributions as, e.g.vertical profiles.Especially in the reconstruction of a 2-D distribution, the error due to interference depends on the constraint and is important, as it describes how independent the different quantities are.This error can be calculated analogously to the smoothing error using the block matrix of both concerning profiles (Sussmann and Borsdorff, 2007) if a covariance matrix of each interference species is known.But more work including model studies are necessary for generating of the necessary covariance matrix.
As the dimension of the solution vector x in Eq. ( 6) is three times larger than the measurement vector y, it is possible to obtain the same fit either by adjusting only one component of the wind velocity (e.g.x or y) or the source flux.The choice of the constraint distributes the information to the different quantities.Therefore, the different retrieved quantities are dependent not only on their own constraint, but also on the constraint of the other quantities.The Optimal Estimation Theory requires the a priori covariance matrices of all of these quantities and distributes the information automatically.In this work the constraint is adjusted empirically and we can just argue that the dependence between the source flux and wind velocity is minor, as the results of both shown in Figs. 3 and 4b have very different structures.However, if the constraint of the wind velocity field does not allow for variations on smaller scales, an under-constrained source flux field might compensate this so that the missing small structure appears now in the source flux field and the simulation fits the observation well again.This would diagnose a direct dependence between both the source flux field and the 2-D wind velocity field.However, the retrieval of the source flux field is plausible.

Total error
The total error in the mean emission rate is estimated from the independent errors in (a) the possible systematic error of the slant column retrieval (≈ 20 %), (b) the geometrical error (≈ 21 %) (Sect.4.  Stremme et al., 2012a), an independent wind field retrieval can be performed and compared, as shown in Figs.8a and b.
The retrieved wind speeds perpendicular to the line-of-sight from SiF 4 and SO 2 were found to be 0.8 m s −1 and 0.7 m s −1 , respectively.This fits to 0.7 ± 0.3 m s −1 , the perpendicular wind component towards east of the average and standard deviation of the wind between 5700 and 5900 m a.s.l., given by the radiosonde launched at 06:00 LT at the Mexico City airport.In the image the plume centre is seen at an elevation angle around 5 • , so that the estimated height in the distance of 12 km is around 5800 m a.s.l.
As the SiF 4 column retrievals have large random errors, the velocity retrieval is rather strong constrained and only the average velocity was determined.However, the form of the SiF 4 plume and the propagation determined are confirmed using the SO 2 plume as tracer, indicating that the volcanic gases first rise and then drift to the east in this particular example (left from the image).The horizontal velocity is rather small, and the initial ascent of the gas is visible.Unfortunately, the wind direction close to the source is not perpendicular since the velocities are strongly constrained towards a mean average velocity.

SO 2 emission rate
The estimated emission rate on the night of 1 to 2 December 2007 (22:43 LT), shown in Figs. 3 and 5, is found to be around 1000 t d −1 (with a maximal value of 1200 t day −1 ).The result of this "snapshot" is a typical value for the Popocatépetl and only slightly smaller than the average emission rates found in the literature.SO 2 emission of Popocatépetl has been estimated for different times and periods.An average emission rate of 2000-3000 t day −1 is reported by Delgado-Granados et al. (2001) for the preeruptive phase, and an average with standard deviation of (2450 ± 1390) t of daily SO 2 emission was found by Grutter et al. (2008a) during an intensive measurement campaign in March of 2006.The average of the daily SO 2 emission during a two-year period (2005)(2006)(2007) was calculated by the NOVAC project and found to be 1503 t day −1 with a standard deviation of 1154 t day −1 (Rivera, 2009).
The reconstructed wind velocity can be compared with the velocity measured by the radiosondes launched at 06:00 LT near the airport.The radiosonde reports a wind speed of 2.1 m s −1 and a wind direction of 69 • N in the altitude 5450 m.a.s.l.As Altzomoni is located at 343 • to the north of the crater, an exact perpendicular wind direction would be 73 • , which is almost given by the radiosonde measurements.This is a necessary condition to compare wind speeds directly, otherwise the wind velocity has to be projected to the perpendicular plane.The mean wind speed retrieved with this experiment (2.0 m s −1 ) is consistent with the radiosonde measurement.
This work presents a method for estimating the propagation of a volcanic plume from measured thermal IR gas columns by constructing a 2-dimensional wind field.The reconstruction is based on the equation of continuity and uses a pair of sequential 2-D distributions of the column densities.The column densities are provided by a scanning infrared gas imaging spectrometer, as described in Part 1, but the algorithm can be applied to other imaging devices as well.The wind vectors can be displayed online over the gas-distribution sequences as an animation (see Supplement), and the overall propagation velocity serves to calculate emission rates from the volcano.
The reconstruction method uses optimal estimation and a Tihkonov-like regularisation, allowing for the definition, calculation and the use of averaging kernels in analogy to atmospheric profile retrievals.The algorithm has been tested on volcanic plumes measured on several occasions during the 2007-2009 period from Popocatépetl in Central Mexico.Plume propagation velocities retrieved independently from SiF 4 and SO 2 were presented in one example and found to be self consistent.
It is shown that with favourable conditions, the average velocity from SO 2 column densities could be retrieved with this algorithm and confirmed by an independent crosscorrelation method and from radiosonde data.In one example the calculated SO 2 emission rate of 1000 t day −1 is a reliable result and consistent with the current knowledge of Popocatépetl's emission rates and variability.This methodology complements other techniques like the mobile and scanning DOAS instruments, but adds the possibility of monitoring also during the night and with an entirely different approach.A longer time series would be desired in order to allow for cross-validation of the calculated emission rates and compare with other methods.Two major improvements regarding the forward model should be addressed in future work: (1) there is an inconsistency because the scanning device is continuously monitoring until the image is constructed, while the forward model treats the images as snapshots, and (2) the 3rd dimension is missing.In principle the forward model could be completed with the missing 3rd dimensions by using (a) the sensitive dependence of the slant column retrievals on the plume distance to the observer and (b) the simple geometrical consideration that a plume moving towards the observer increases.1 is of block-shape:

A.
Krueger et al.: Thermal emission spectroscopy of volcanic gases -Part 2

Fig. 1 .
Fig. 1.A sequence of two consecutive SO 2 column images taken 3 min apart by thermal emission spectroscopy and a continuously scanning device on 17 March 2006.The column amount of SO 2 is given by the false-colour image, warm colours meaning high SO 2 , cold colours meaning low SO 2 column densities.The dashed lines depict the origins of two puffs within the plume detected in frame (a) which has moved to the solid lines in frame, (b).The average wind vector is indicated by the arrow.The wind vectors describing the plume propagation projected in the plane of the image are sought.

Fig. 2 .
Fig. 2.False-colour plot showing the difference of (a) the observed column distribution between two consecutive frames, representing the measured input vector for the reconstruction of the wind field.(b) is the same as (a) but calculated using the forward model K and the solution vector x, consisting of the retrieved vector field and the retrieved source field.The residual (c) is the difference between the plots above (observed -simulated).Due to the choice of S −1 e (see text, Sect.3.4), the residual of the column differences are larger outside of the plume.The location of the Popocatépetl volcano is indicated.

Fig. 4 .
Fig. 4. (a)Diagonal of the regulation matrix concerning the strength of the sources.As the units are arbitrary, they have to be compared with the constraint (R) concerning the wind components and with the weighting of the measurement (S −1 e ).(b) Retrieved sources are also part of the solution vector.As the wind field is constrained, this additional fit parameter improves the result significantly.The clean atmosphere is constrained with a smoothing constraint, while the pixels of the borders and the volcanic body are constrained with an optimal-estimation-like regularisation (R-matrix is diagonal).In this example, 11.7 independent pieces of information are included in the vector, which has 714 elements to describe the sources.As expected, the retrieval finds the strongest sources above the crater.
for the weighting of the measurement taken in each pixel.The regularisation matrix Three strategies might solve a mathematically ill-posed problem: (i) reduction of the resolution of the solution www.atmos-meas-tech.net/6/47/2013/Atmos.Meas.Tech., 6, 47-61, 2013 A. Krueger et al.: Thermal emission spectroscopy of volcanic gases -Part 2

Fig. 5 .
Fig.5.Calculated SO 2 emission (t day −1 ) at different distances from the crater represented by the cross-sections drawn in Fig.3and using the reconstructed wind field.The distance is shown as the time of its propagation after the starting point of the linear trajectory (Fig.3).The solid trace shows the flux calulated using the first image of SO 2 columns and the dashed line is the flux at the same time, but based on the next SO 2 image.The reconstruction of the wind field is done using a constraint which limits its spatial resolution.Especially near the crater, the dispersion of the SO 2 plume is not very well represented by the retrieved wind field and the SO 2 emission is underestimated (more or less below 6 min).A similar effect occurs at the other end of the image.In the centre of the plume, where the linear trajectory is well represented, both SO 2 emission estimations are consistent and oscillate around 1000 t day −1 in this case.

Fig. 7 .
Fig. 7. Diagonal of the averaging kernel of the (a) x-(A x ) and (b) y-wind velocity components (A y ), as well as of (c) the sources (A src ).(I.)The first two retrievals with stronger Tikhonov and with negligible OET (diagonal) constraint, (II.)Final retrieval with negligible Thikonov constraint but stronger OET (diagonal) constraint.As the main wind direction is along the x-axis, A x represents the sensitivity to a deviation from the wind speed (|v| − |v a |).As the main wind direction is along the x-axis, A y represents the sensitivity to a deviation from the wind direction (dv y = |v a |dφ).The sensitivity to sources depends on the constraint (R) and the weighting of the measurements (S −1 e ).
et al.: Thermal emission spectroscopy of volcanic gases -Part 2

Fig. 8 .
Fig.8.Wind fields reconstructed from consecutive measurements of (a) SiF 4 and (b) SO 2 .Even when the retrieved SiF 4 columns have rather large random errors and low precision, an average velocity can be obtained independently from the SO 2 plume retrieval.The retrieved wind speeds perpendicular to the line-of-sight from SiF 4 and SO 2 were found to be 0.8 m s −1 and 0.7 m s −1 , respectively.

A
3.1 Stremme et al., 2012a, Part 1) and (c) the measurement error (≈ 20 %, Sect.4.3.2) of the resulting mean wind speed and emission rate.The overall error of the calculated emission rate, still missing the contribution of the smoothing and interference errors, is around 35 %. wind field retrieval has already been shown in Fig. 3 from SO 2 columns measured during passive degassing periods of the volcano.When another gas can be detected, like SiF 4 during the eruptive event of 18 November 2008 (see Part 1,

Table 1 .
Overview of the retrieval schema: WD -wind direction, WS -wind speed, WV -wind velocity, CC -cross correlation, Aaveraging kernel.lag describes the time lag resulting from the CC, t frame is the time difference between two images. t

Krueger et al.: Thermal emission spectroscopy of volcanic gases -Part 2
Abbrev.description j former ij trace gas flux in pixel ij at the time when the first image is taken i latter ij trace gas flux in pixel ij at the time when the second image is taken E former A.