Interactive comment on “ Monitoring volcanic ash cloud top height through simultaneous retrieval of optical data from polar orbiting and geostationary satellites ” by K

Volcanic ash cloud-top height (ACTH) can be monitored on the global level using satellite remote sensing. Here we propose a photogrammetric method based on the parallax between data retrieved from geostationary and polar orbiting satellites to overcome some limitations of the existing methods of ACTH retrieval. SEVIRI HRV band and MODIS band 1 are a good choice because of their high resolution. The procedure works well if the data from both satellites are retrieved nearly simultaneously. MODIS does not retrieve the data at exactly the same time as SEVIRI. To compensate for advection we use two sequential SEVIRI images (one before and one after the MODIS retrieval) and interpolate the cloud position from SEVIRI data to the time of MODIS retrieval. The proposed method was tested for the case of the Eyjafjallajokull eruption in April 2010. The parallax between MODIS and SEVIRI data can reach 30 km, which implies an ACTH of approximately 12 km at the beginning of the eruption. At the end of April eruption an ACTH of 3–4 km is observed. The accuracy of ACTH was estimated to be 0.6 km.


Introduction
Volcanic ash cloud top height (ACTH) is an important parameter in aviation route planning (Sparks, 1997) as well as in various volcanological models for estimating a volcano's mass eruption rate (Mastin et al., 2009;Stohl et al., 2011).ACTH can be estimated using pilot reports, ground-based, airborne, and spaceborne measurements.Ash clouds have been observed from the ground by weather radar in a number of instances, e.g. Rose et al. (1995) and Lacasse et al. (2004).During the last couple of years, observations of volcanic ash by ground lidar, e.g.Gasteiger et al. (2011) or Hervo et al. (2012), received increasing attention.Although these are excellent tools, all ground-based measurements are, however, restricted to their exceptional spatial and temporal availability.
Compared to ground-or airborne-based observations, satellite remote sensing of ACTH has the great advantage of being available all the time.Geostationary satellites even provide a temporal resolution of some minutes.The most reliable data are retrieved using active instruments like spaceborne lidar (e.g.CALIOP, see Sect. 2).However, only one such instrument is operating nowadays, making it unsuitable for operational monitoring of volcanic plumes because of the long revisit time.A more common method is based on the measurements of cloud brightness temperature, cloud emissivity, and atmospheric temperature profile.Because of the difficulties obtaining all parameters accurately, this method can lead to considerable ACTH errors.Improvements to this method can be made using radiative transfer models, as in e.g.Stohl et al. (2011).As another alternative we propose here a new solution for ACTH retrieval based on a photogrammetric method.
In the following, we first survey available spaceborne measurements, including a detailed review of photoclinometric and stereoscopic ACTH measurements.The proposed photogrammetric method is described next (Sect.3) and case study results are analysed thereafter (Sect.4).We finish with a discussion of the method's accuracy in Sect. 5.

ACTH based on active sensors measurements
The Cloud-Aerosol Lidar with Orthogonal Polarization (CALIOP) flies aboard the CALIPSO satellite (NASA, 2011).This instrument has a horizontal resolution of 333 m and vertical resolution between 30 and 180 m, depending on the distance to the ground.It has a revisit time of 16 days and a swath width of only one kilometre and performs only nadir measurements.CALIOP has already been used for volcanic ash cloud monitoring during, e.g.eruptions of Chaiten 2008 (Carn et al., 2009), Kasatochi 2008 (Karagulian et al., 2010) and Eyjafjallajökull 2010 (Stohl et al., 2011).
Another active instrument that might be used to determine ACTH height estimates is the Cloud Profiling Radar (CPR) aboard CloudSat (Colorado State University, 2011).It has lower vertical (240 m) and horizontal resolution (1700 m) than CALIOP, a revisit time of 16 days, and just like CALIOP provides only nadir measurements.CloudSat and CALIPSO fly in formation (same orbit with 15-s time delay), thus their results can be easily combined.
Because of the orbit characteristics the CALIOP and CPR instruments are prone to miss a volcanic eruption.Thus, they are not appropriate for operational monitoring, but are very useful for validation and supplementary information in case data are available.

ACTH based on brightness temperature
Satellite brightness temperature (BT) methods have been used several times to estimate ACTH (Oppenheimer, 1998;Prata and Grant, 2001;Tupper et al., 2004).This method compares BT retrieved from the ash cloud (normally utilizing the 11 µm window channel) against the local vertical atmospheric temperature profile.The height at which the retrieved BT matches the atmospheric temperature profile is considered to be ACTH.
Several factors limit this technique (Oppenheimer, 1998).The first limitation is the assumption that the ash cloud top is in thermal equilibrium with the ambient air.If the ash cloud overshoots its thermally equilibrated level, or if the ash cloud still has sufficient energy to develop to higher altitudes (after the time of the satellite image), the BT method does not work well.
The second limitation is the assumption that the ash cloud emissivity equals 1.The actual emissivity is, however, not known.BT of thick ash clouds will under this assumption closely approximate the true BT.But in the case of a dilute cloud the spaceborne instrument will detect radiation from beneath the ash cloud, effectively lowering the heights.
Third, BT methods also need an accurate atmospheric temperature profile.This profile depends mostly on season and latitude, but even smaller variations can lead to significant ACTH errors.Additional inaccuracy is caused by the instability of the temperature profile near the tropopause (Prata and Grant, 2001).In case of the Eyjafjallajökull eruption, the tropopause is not higher than 10 km, thus the ACTH of 9 km estimated from BT is already questionable.In tropical regions there are additional problems due to the so-called double tropopause (Randel et al., 2007).

ACTH based on absorption by trace gasses
CO 2 absorption technique is a common method for estimating the cloud top height (Chang et al., 2010).It usually requires measurements at 11 µm as well as at 13-15 µm.The CO 2 absorption technique makes use of the atmosphere becoming more opaque due to CO 2 absorption as the wavelength increases.Therefore, the radiances measured in spectral bands between 11 and 15 µm are sensitive to a different layer of the atmosphere.To estimate the cloud top height, the differences between cloudy and cloud-free areas are rationed in a CO 2 band and a thermal infrared window band.These values are then compared with the theoretical values to find the optimal matching height.The method yields the best accuracy for the upper troposphere.The problem is that the required bands are not available on many operational instruments (e.g. the Advanced Very High Resolution Radiometer (AVHRR), or even on the recently launched Visible Infrared Imager Radiometer Suite -VIIRS).The method has already been tested for ACTH estimation (Richards et al., 2006); isothermal regions and clouds at or near the tropopause cannot be retrieved using this algorithm because of the instability of the temperature profile near the tropopause.Another problem is the underestimation of ACTH because the method retrieves a radiative height ("height of the cloud centre").
Height of (aerosol) clouds can be estimated from reflectance ratio measurements in the O 2 absorption A-band at 0.76 µm (Dubuisson et al., 2009;Preusker et al., 2007).The ratio is computed from a reflectance of a spectral band, strongly attenuated by O 2 absorption, and from the reflectance in a second spectral band with minimal attenuation.For a given surface reflectance, simple relations have been established between the reflectance ratio and the altitude of an aerosol layer as a function of atmospheric conditions and the geometry of observation.The method has been tested for POLDER (aboard ADEOS satellite) and MERIS (aboard ENVISAT satellite).The simulations show that the method is only accurate over dark surfaces when aerosol optical thickness is larger than 0.3.The same methodology is currently being tested also for ACTH estimation (ESA, 2011).

ACTH based on backward trajectories modelling
ACTH based on backward trajectories modelling correlates cloud movement with atmospheric winds (Eckhardt et al., 2008;Oppenheimer, 1998;Tupper et al., 2004).This method takes advantage of vertical wind profiles measurements -the horizontal wind component at any given altitude is usually unique in its combination of direction and speed.Therefore, airborne ash moves in this direction with the speed of the prevailing wind.If the direction and speed of the airborne ash is determined from the satellite measurements, an estimate of ACTH is possible by matching the cloud movement with the corresponding wind movement.The method usually uses data from low spatial resolution satellite instruments with revisit times of about 12 h, like for instance the Infrared Atmospheric Sounding Interferometer -IASI (Klüser et al., 2013).In general, any kind of instrument can be used as long as it provides the data required by the model.It is the model that "moves" the detected ash cloud back to its source, thus the quality of the model and meteorological data (wind velocity and direction as a function of height) are of highest importance.

ACTH based on shadow lengths and photoclinometry
A trigonometric way of estimating ACTH is to measure the ash cloud shadow lengths (Glaze et al., 1989;Prata and Grant, 2001) on daytime images.The processing is more difficult if the method is used over rough terrain as in this case the terrain elevation has to be considered.More importantly, the value of the shadow lengths usually corresponds to the margins of clouds, rather than their highest central points.To measure the top of the cloud height, the sun should be close to the horizon.But a large zenith angle has a consequence of low contrast, thus the shade mask is unreliable.A possible improvement is the combination of shadow length and temperature determinations using a photoclinometric study (Glaze et al., 1999;O'Hara and Barnes, 2012) which showed that the shadow technique alone systematically underestimates the ash cloud height by up to 30 %.

ACTH based on stereoscopy
Parallax shifts of meteorological clouds or volcanic ash clouds are not obvious to many satellite data users because they mostly use imagery from a single satellite.However, when data from two different satellites are compared, the parallax shift becomes obvious -it is larger for higher clouds and plumes.Although proposed already before (Ondrejka and Conover, 1966), it was Hasler (1981) and Hasler et al. (1983) who drew attention to this by using two geostationary satellites (GOES) for stereoscopic measurements of cloud-top heights.Stereoscopy from geostationary satellites was also used for monitoring cumulonimbus clouds (Hasler et al., 1991), for research on convective clouds (Fujita, 1982;Mack et al., 1983), and for developing climatologic cloud monitoring systems (Wylie and Menzel, 1989).Because of the malfunction on GOES 6, the development in this field stopped until GOES 8 and GOES 9 were launched (Wylie et al., 1998).
In the studies listed above the cloud top was determined from a pair of GOES satellites.Seiz et al. (2007) used Meteosat images to retrieve cloud heights.The Meteosat-5/-8 HRV and Meteosat-5/-7 combination were tested.The resulting accuracy is about 1000 m which is worse than the estimated accuracy by GOES of 500 m (Hasler et al., 1983).The combination of Meteosat-5/-8 TIR data was used also for the eruption of Karthala (November 2005); Carboni et al. (2008) estimated the ACTH to be between 11 and 15 km above sea level.
An alternative to a pair of geostationary satellites are instruments with multi-angle observation possibilities.Prata and Turner (1997) proposed an algorithm for determining cloud heights based on the ATSR (aboard ERS satellites) forward and nadir views of clouds.The method has an accuracy of approximately 1000 m.It was used to determine ACTH for the Mt.Ruapehu eruption in 1996.Data were retrieved from ATSR-2.The heights on 8 July were about 8000 m for downwind points, and 5500 m for points near to the vent.Muller et al. (2007) showed that ATSR-2 cloud-top heights depend on the wavelengths used in the stereoscopy.They proposed that a combination of visible and thermal bands could yield information on multi-layer clouds.In addition, Along Track Scanning Radiometer (ATSR) was used also for cloud detection in polar regions (Cawkwell et al., 2001).
In 1999 the Terra satellite was launched.Aboard are, among other instruments the instruments, the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) and Multi-angle Imaging SpectroRadiometer (MISR).Both are capable of stereoscopic height measurements.Comparing different instruments it was shown that ASTER stereo cloud top heights are on average 1000 m higher than Moderate Resolution Imaging Spectroradiometer (MODIS) BT heights (Genkova et al., 2007).Unfortunately, ASTER has only a 60-km-wide swath and thus a long revisit time.
In the case of MISR, the revisit time is shorter (9 days, swath width of 360 km), thus it is more appropriate for use in volcanology.MISR was utilized to retrieve ACTH, optical depth, type, and shape of the finest particles of two highly explosive eruptions occurring on Mount Etna in 2001 and 2002 (Scollo et al., 2010(Scollo et al., , 2012)).For the Eyjafjallajökull eruption in 2010, it was also used to validate model estimated ACTH (Stohl et al., 2011) and plume particle-type characteristics (Kahn and Limbacher, 2012).
An interesting study was the use of the infrared spectral imaging radiometer that flew as part of mission STS-85 on the space shuttle Columbia in 1997 (Lancaster et al., 2003).A method for computing cloud-top height with a precision of 620 m based on multispectral stereo measurements was developed.The study is interesting because the results are compared with coincident direct laser ranging measurements from the shuttle laser altimeter (the laser altimeter mean heights were about 100 m lower).Another study using the same datasets (Manizade et al., 2006) showed that the accuracy of the height can be improved if the data are first segmented according to brightness temperature.
The combination of two different instruments, one aboard a satellite in geostationary and the other aboard a polar orbiter has also been suggested, but never really exploited (Hasler et al., 1983).With several new instruments being in orbit since 1983 this combination offers new exciting possibilities that are detailed in the following section.

Proposed method of ACTH based on a parallax
The proposed method estimates ACTH based on the parallax of the data retrieved nearly simultaneously by two instruments.Figure 1 shows the principle of the method: satellites in the polar and geostationary orbit do not observe an ash cloud at the same viewing angle.This causes parallax when comparing both images.The parallax depends on the differences in the viewing geometry and the height of the ash cloud.As the viewing geometry is defined by the coordinates of the retrieved data and the coordinates of both satellites it is easy to retrieve ACTH.
In our case study (Eyjafjallajökull eruption in April 2010; see Sect. 4) we use a combination of MODIS aboard the Terra and Aqua satellites (polar orbit) and the Spinning Enhanced Visible and InfraRed Imager (SEVIRI) aboard the Meteosat Second Generation 2 (MSG2) satellite in a geostationary orbit.Both instruments are multispectral, but we limit the scope of this study to data from the visible spectrum, allowing for ACTH estimations only during the day.The reason for this limitation is the spatial resolution of the retrieved data.The resolution has (besides parallax) the largest influence on the photogrammetrically estimated ACTH.The HRV band of SEVIRI has a nominal resolution of 1000 m (in nadir) and MODIS bands 1 and 2 have nominal resolution of 250 m (in nadir; we used only band 1).These two resolutions are significantly higher than the resolutions of the accompanying thermal bands (3000 m for SEVIRI and 1000 m for MODIS).
The proposed method of ACTH estimation consists of three main steps.In the first step we aggregate MODIS data to the SEVIRI spatial grid.The second step is automatic image matching.In the third step, lines of sight connecting observed points of both satellites are generated; the intersection points of SEVIRI and MODIS lines of sight are then used to estimate ACTH.

Data pre-processing
To be able to perform automatic image matching, it is necessary to pre-process data so that MODIS and SEVIRI datasets are comparable.In the previous retrievals of meteorological cloud-top height (Hasler, 1981), both images from GOES were projected to a standard map projection.All such transformations usually decrease the accuracy of the results.For the case study presented in Sect.4, SEVIRI, with its poor spatial resolution, is the parameter that limits the accuracy of the results, thus we decided to leave SEVIRI data in their own grid system.MODIS data have much better spatial resolution, thus they can be projected to the SEVIRI grid system without significantly influencing the resulting accuracy.The SEVIRI grid system is a normal geostationary projection, known also as near-side general vertical perspective projection.Therefore, the easiest way to transform a MODIS image to the SEVIRI grid is to use the just mentioned map projection.The MSG2 satellite's position is, however, not stable.The satellite might change its position even for a hundred kilometres.Such changes in the satellite's position might lead into a positional error on the order of one pixel.Therefore, we prefer to not use a map projection to generate a MODIS dataset in the SEVIRI grid system, but an alternative procedure (example in Fig. 2): -Estimate which MODIS pixels contribute to each SE-VIRI pixels.
-Estimate the distance of each MODIS pixel from the corresponding SEVIRI pixel centre.
-Convert these distances to weights considering the SE-VIRI point spread function (Deneke and Roebling, 2010).
-Normalize the weights in such a way that their sum equals one for each SEVIRI pixel.
-To get the final MODIS based dataset in the SEVIRI grid, multiply all MODIS pixels with their weights and sum these values for each pixel of SEVIRI's grid.

Image matching
The goal of the satellite image matching step (called also image-to-image cross-correlation) is to accurately identify point pairs between two satellite images.This might be difficult if the images are not retrieved by the same instrument.
The problem involves different resolutions, different viewing geometries, and different instruments' response functions.In addition, the appearance of the same object in two different images might contain a large illumination variation, and thus the local descriptors of the same feature point are different.A number of automatic image matching approaches have been proposed to solve these issues.The methods can be categorized as area-based or feature-based methods (Evans, 2009;Goncalves et al., 2011;Hasan et al., 2010;Huo et al., 2012;Teke and Temizel, 2010).Scale invariant feature transform (SIFT) (Goncalves et al., 2011) and speeded-up robust features (SURF) (Evans, 2009) are becoming popular feature-based methods for camera image registration due to scale and rotation invariance of the detector and the distinctiveness of the descriptor.Simple matching based on Euclidean difference between SURF descriptors is not sufficient.Several authors have proposed geometrical constraints such as scale and orientation restriction (Teke and Temizel, 2010) to improve the number of correct matches for remote sensing images.Our experiments proved that this was not enough in order to obtain a sufficient number of correct matches.We thus applied a more strict location relationship among neighbouring feature points, which allows only for small local position errors within a global affine transformation.Rather than using a single global mapping function for registration, we use local image patches, which are registered using an affine model.Further experiments showed that it is possible to find matching features on the volcanic ash cloud.The problem is that only a small number of points are selected for further processing.In addition, it is difficult to control, which features are really chosen -we cannot be sure that the chosen features represent the highest parts of the volcanic ash cloud.
Although the feature-based image matching (SURF) performs good (results correspond well to results of manual points selection), we also tested the area-based image matching method.Such an approach allows the detection of more pixel pairs than the SURF-based method.Many other approaches were used in previous studies (Hasler et al., 1991;Löfdahl, 2010;Ramapriyan et al., 1986;Scambos et al., 1992;Wu et al., 2012), some of which are based on the frequency domain, and some on the image domain.As pointed out by Prata and Turner (1997), the frequency domain methods based on the fast Fourier transform makes sense only when working with large images with a constant shift all over the image.Within this study we were not comparing whole images, but just parts of it using a moving window analysis.This was necessary as the case study area (see Sect. 4) is large; in such a large area the clouds do not necessary behave the same in all parts of the image.Thus, it is necessary to determine local shifts by a moving window analysis.We used the same procedure as already described by e.g.Scambos et al. (1992) or Prata and Turner (1997).The method contains the following steps: -Based on the pixel being currently processed (with image coordinates C, L) select a reference subset (with a number of columns n C1 and a number of lines n L1 ) of the first image.This is shown in Fig. 3 (left): pixel C, L is in red, the reference frame is yellow.
-In the second image select a search subset (with a number of columns n C2 > n C1 and a number of lines n L2 > n L1 ) that is centred over the pixel having the same image coordinates (C, L).This step merely reduces the possible extent of a data subset from both images to be compared.Thus, it is theoretically not necessary but it reduces the computing time.The size of this search subset should be large enough to detect the largest expected shifts.Figure 3 (right) shows the search subset (light orange) in the second image centred at the pixel C, L.
-Within the search subset define a moving window (with a number of columns n C1 and a number of lines n L1as large as the reference subset in the first image).The moving window is shown in yellow in Fig. 3 (right).
-Move this window across the search subset, and for each possible shift of the moving window compute the  (1) where CI is the correlation index between both subsets, DNm i,j and DNr i,j are digital numbers of the moving window (second image) and reference subset (first image), µm and µr are the mean values of reference subsets and the moving window set, i and j are shifts between the central pixels of the reference subset and the moving window.
-For the central position C, L of the reference subset, return the shift (the position of the moving window centre) where computed CI is the largest.In Fig. 3 (right), the centre of the moving window with highest CI is red, and the centre of the corresponding reference frame in the first image in dark orange.Matching shift between them is noted as dC and dL -the cloud is in this example shifted two columns to the right and three lines to the bottom.If n C and n L are odd numbers, then the maximal possible shifts are dC = (n C2 − n C1 )/2 and dL = (n L2 − n L1 )/2.
The use of image matching is restricted to an estimation of image shifts between a pair of images.The shifts can be used to calculate those image coordinates in the second image that correspond to the points of the first image.In the proposed procedure of the ACTH estimate, three images are always used: one MODIS, a SEVIRI image acquired just before, and a SEVIRI image acquired just after the MODIS image was retrieved.Therefore, image matching has to run twice to find matching points in all three images.
Results of image matching depend on the size of the search area and moving window.A large moving window can detect large features, but it usually fails to detect small features.In contrast, a small moving window detects small features but generates a lot of noise in the results.The appropriate optimization is image matching over image pyramids.We consider image pyramids as a multi-resolution representation of the original image (Anderson et al., 1984).Each higher pyramid is merely a regridded lower pyramid.The lowest pyramid is the original image.With each higher pyramid the resolution decreases for a chosen pyramid factor -in our case a factor of 3. Thus, the data of a higher pyramid are always averaged within a 3 × 3 pixels large area of its lower pyramid.In our case we use three pyramids.For instance, if we had original data of 270 by 540 pixels, we would aggregate them first to a pyramid of 90 by 180 pixels and in the following step to a pyramid of 30 by 60 pixels.
Shifts are first estimated over the image pyramid having the coarsest resolution.A search subset is defined (Fig. 4, point a); the central points of the search subset coincide with the central points of the reference subset.Image matching estimates the shift using moving window analysis (Fig. 4b and  c).These shifts are then refined on a lower pyramid.First, they are multiplied by the pyramid factor and are regridded to the resolution of a lower pyramid (Fig. 4, point d).These shifts then define new central points of each search subset in this pyramid.This procedure is repeated until the original data (the finest image pyramid) are processed.Such a scheme makes results more reliable and the computation is significantly faster.If CI was lower than 0.7 between the compared datasets in a coarse pyramid, the estimated shifts were considered unreliable and set to 0. For all three levels of image pyramids the moving window was 7 × 7 pixels large and the search area 13 × 13 pixels large.Considering the geometry in the case study area, the size of the search subset and the chosen three image pyramids restrict the maximum ACTH to about 30 km.

Intersection of lines of sight for ACTH estimation
After the triple of the corresponding image coordinates (from MODIS and both SEVIRI images) is known, it is easy to compute their coordinates in a global coordinate system.The data required for the estimate of the longitude and latitude is usually part of the satellite dataset or can be obtained from their metadata.Also, the satellite's ephemerides are usually a part of the metadata, which allows computing the satellite's position at the time of the data retrieval.Having these data, it would be possible to estimate ACTH following Prata and Turner (1997) from appropriate zenith angles.Here we propose another solution based on vector algebra.
-The first step is conversion of all available coordinates into a geocentric Cartesian coordinate system (see Appendix A).
-The effect of possible advection of the eruption cloud between the MODIS and the SEVIRI images is considered for each pixel triple: the coordinates of a virtual SEVIRI pixel are interpolated from position of both SE-VIRI pixels to the time of MODIS retrieval.
-We define parametric equations of 3-D lines connecting coordinates of the virtual SEVIRI pixels with the position of the MSG2 satellite ("SEVIRI lines") and corresponding lines connecting coordinates of the MODIS pixels with the position of the Terra/Aqua satellite ("MODIS lines").
-The solution of the following linear system gives the intersection of such a line pair: where [x,y,z] M and [x,y,z] S are the positions of the MODIS aboard Terra/Aqua and SEVIRI aboard MSG2., [v x ,v y ,v z ] M and [v x ,v y ,v z ] S are direction vectors of MODIS and SEVIRI lines, t M and t S are the unknowns defining the point of intersection.
-The system in Eq. ( 2) is overdetermined, thus it can be solved by a least-square technique.The geocentric www.atmos-chem-phys.net/13/2589/2013/Cartesian coordinates of the intersection are then converted back to the geographic coordinate system: longitude, latitude, height above ellipsoid -i.e.ACTH (see Appendix A).
MODIS and SEVIRI lines never intersect because the data are not continuous but discrete pixels.The lines rather pass each other.Thus, the Eq. ( 2) does not search for the real intersection but for the pair of closest points on the corresponding lines.ACTH can then be estimated from one of these two points or as their average.The advantage of this ACTH esti-mate compared to the solution from zenith angles (Prata and Turner, 1997) is the possibility to estimate the intersection quality.It can be described by the distance between MODIS and SEVIRI lines; if it is small, the accuracy of ACTH is high.In the following, we use the expression "intersection distance" to express this value.

Case study: eruption of Eyjafjallajökull in April 2010
An explosive eruption under the Eyjafjallajökull glacier on Iceland started on 14 April 2010.The interaction of magma and melted water amplified the explosive activity, which resulted in emissions of fine ash and volcanic gas into the atmosphere (Gudmundsson et al., 2010).Intense eruptive activity continued for five days.During this period, prevailing meteorological conditions resulted in ash transport directly towards Europe (Fig. 5).The eruption strength increased again in May 2010.
In this case study we focus on the eruption period 15-19 April.The eruption featured an initially phreato-magmatic phase followed by a magmatic phase (Gudmundsson et al., 2010).However, volcanic eruptions undergo different phases of activity even during short periods of time on the order of seconds to minutes (Scharff et al., 2012).This can be seen in Fig. 6, which presents brightness difference temperature (BTD) from SEVIRI bands 9 and 10, with negative BTD being an indicator of ash presence in the atmosphere (Prata, 1989).The reddish-coloured plume over the North Sea spreading towards Denmark is airborne ash.On its western end it is covered by a dark-blue-coloured cloud that is not recognized as airborne ash.BTD for this darkblue-coloured cloud exceed +10 K, being indicative of a meteorological cloud.However, analysis of the animated data shows that this blue-coloured plume is also of volcanic origin.It spreads from Eyjafjallajökull across the North Atlantic and the North Sea and then it turns to Northern France.During the phreato-magmatic eruption, more water is injected into the atmosphere compared to a regular eruption.As the eruption cloud rises and cools, the water starts to condensate and eventually will form ice.In this case, condensation and ice formation are certainly enhanced by ash, which serves as nuclei for water and ice condensation.Therefore, one can as-sume from BTD that the blue plume in Fig. 6 consists not only of ice, but of ice with ash.This blue plume is located higher than the plume in red.Because of this height differences the dispersion of both plumes is different.The red part moves first to the east and eventually to the south-east and the blue part to the south-west.
Figure 7 shows a SEVIRI image taken on 17 April at 13:12 UTC (Fig. 7a) that was retrieved some minutes before the MODIS image.South of 60 • N a large area of meteorological clouds covers the area of the North Atlantic, Northern Ireland and Scotland.These clouds reach heights from 4 km up to over 9 km.After the initial part of the eruption they moved over the North Atlantic, covering Iceland on the evening of 15 April.This restricted monitoring of volcanic ash over Iceland by thermal remote sensing.Once the clouds went away over Iceland in the beginning of 17 April, a small ash cloud could be observed again.This can be seen in Fig. 7a stretching from Iceland first to the south and then turning to the east.Note that the majority of the ash cloud is positioned north of 60 • N. Furthermore, a very high cloud can be observed over the North Sea.From the SEVIRI image (Fig. 7a) the ash cloud does not seem thick compared to the previously discussed meteorological clouds.Its origin is difficult to determine.BTD is highly positive, indicating a meteorological cloud.However, analysis of the animation shows that the cloud first appeared above Scotland.As it remained moving in the similar direction as the first ash cloud reaching Scotland, we assume that this is also an ice cloud originating from the volcano that was for some time "hidden" by higher clouds.
Figure 7b shows the cross-correlation index (CI; Eq. 1) between the shown SEVIRI image (Fig. 7a) and corresponding MODIS image retrieved at 13:15 UTC.As expected it is very high along the coastal areas having large contrast between the land and the sea.The lowest CI is within large land areas and the open sea not covered by clouds.Similar characteristics can be observed in cloudy areas, where clouds contain some edges and texture CI is high as well.CI is also well correlated with estimated column shifts (Fig. 7c) and line shifts (Fig. 7d) between both datasets.Where CI is low, sudden changes of shifts are possible.Especially shifts in the line direction are influenced mostly by the parallax between MODIS and SEVIRI data (Fig. 7e).The highest parts of the clouds have a parallax larger than 30 km.This is a shift of MODIS data for about 10 pixels to the north in the case study area.Such a shift can be easily observed even by visual comparison of MODIS and SEVIRI data.
Column shifts are in the case study of lower importance for parallax estimation, but they are highly correlated to the wind velocity.The shifts shown in Fig. 7c-d correspond only to the SEVIRI data retrieved before MODIS retrieval.For an estimate of the wind velocity, the SEVIRI data after MODIS retrieval are needed as well.We do not show them here, but in the case study, the column shifts between both SEVIRI images can reach over 20 pixels.This corresponds to the  velocity of over 80 km/h.The wind velocity can in general influence the estimation of the parallax if the wind direction is the same as the SEVIRI line of sight.If it were ignored in such a case, the cloud advection would be a major source of error.
Figure 7f shows the intersection distance between MODIS and SEVIRI lines of sight (solution of Eq. ( 2) using leastsquare technique); these values can reach distances of some kilometres.Assuming that the accuracy of image matching is half a pixel, then this intersection distance should not exceed this value for a reliable solution.Here we are interested only in the north-south direction because this is the direction in which the parallax mainly extends.The pixel size in the north-south direction is about 3 km in the case study area, thus the intersection distance should here not be larger than 1.5 km.Therefore, many pixels with estimated ACTH (Fig. 7g) are not reliable.It is possible to filter all these pixels (Fig. 7h), but enough pixels remain so that a good image of the cloud topography is possible.The highest pixels represent clouds over the North Sea that reach heights of over 10 km.Another signal shows a meteorological cloud spreading west to east just south of 60 • N. As mentioned before, these are the remains of the cloud that was covering the whole area of Iceland from late 15 April to early 17 April.Its height is mostly between 6 and 9 km.Just north of it, an ash cloud stretching back to Iceland can be observed.It is much lower than the ash cloud was on the first day of eruption but the highest pixels still reach heights over 6 km.
Another example is shown for the first phase of the eruption on 15 April (Fig. 8a).The plume stretches from Iceland over the North Sea north of Scotland towards Norway.The plume can be easily detected directly from the visualized ACTH as it is much higher than the surroundings.The plume ACTH ranges mostly between 5 km and 8 km, but some pixels reach almost 12 km.The plume is, however, difficult to recognize a few days later on 19 April as the eruption intensity has already decreased by that time (Fig. 8b).The ash plume reaches height of over 3 km (orange pixels south of Iceland), which is approximately as high as the surrounding clouds.
Aside from the spatial extent of the cloud it is also interesting to examine the temporal evolution of the ACTH.The analysis covers the ACTH of volcanic ash clouds within the area covering Iceland, Norh Sea, Ireland and UK (the extent shown in Fig. 7). Figure 9 shows how the ACTH vertical distribution changes over the case study period (see Appendix B for the list of used datasets).Figure 9 shows how many ash cloud pixels (in percent) are situated in each height class with 500-m vertical resolution.The maximum ACTH values reached almost 12 km, which corresponds well to the estimation by Stohl et al. (2011).London VAAC reported higher ACTH of over 16 km (Fig. 5), but these values are probably exaggerated due to air traffic safety.
Figure 9 contains also the ACTH estimations from ground radar measurements at Keflavik airport (see Fig. 5 in Marzano et al., 2011).Both datasets show that ACTH tends to get lower over time as the eruption intensity decreases.We note that the radar's vertical resolution is about 2.5 km for the case study geometry.This is a reason for the radar data to oscillate significantly.In addition, the radar data range is limited to 240 km away from the radar position, so most likely not more than 100 km away from the volcano.See Fig. 1 in Marzano et al. (2011) to compare the radar range with our case study area shown in Fig. 7. Our measurements in Fig. 9 contain also data that are even 1000 km away from the volcano.The case study area was not restricted to the same vicinity of the volcano as by radar measurements.This is because too few satellite measurements are then available, making the vertical distribution unreliable.Therefore, the ACTH distribution in Fig. 9 presents not only ACTH at the source, but also during long-range advection.

Discussion
The transport of volcanic ash in the atmosphere and its deposition on land and in the oceans is of significant importance for climate development (Langmann et al., 2010;Robock, 2000).A direct consequence of the airborne ash is usually a closure of the air space for air traffic.Thus, it is very www.atmos-chem-phys.net/13/2589/2013/Fig. 9. Vertical distribution of ACTH for the period between 15 and 19 April.All pixels are classified according to their ACTH in one of the 500 m height classes.Colours represent the relative amount of pixels within a class (in percent) -dark colours show that there are many pixels in a height class.Because of two daytime MODIS overpasses, two ACTH estimations can be processed per day.On 16 April only one estimate is presented because of the unfavourable satellite's swath characteristics.The highest measured heights by a radar located at Keflavik airport (Marzano et al., 2011) are shown for the comparison.Please note that the radar height only reflects the eruption column height near the vent (e.g. at max 100 km away from the vent) and our ACTH cover the whole North Atlantic region.
important to monitor the horizontal and vertical dispersion of volcanic ash during an eruption.The observations of the vertical ash dispersion are still limited to radar and lidar observations, but more methods are available for monitoring the top of ash clouds.One of them is the method presented here.In the following we discuss its accuracy that depends on -ACTH itself, -accuracy of image registration (image positional accuracy), -accuracy of satellite's position, -accuracy of image matching, and -accuracy of the wind velocity.
The influence of ACTH itself on its accuracy can be ignored for satellite remote sensing.The ratio between ACTH and the satellite height above the Earth is too small to significantly influence the intersection between SEVIRI and MODIS lines of sight.The accuracy of the satellite's position and positional accuracy of an image is out of control for a general user of remote sensing data.The position of a satellite is not a significant problem for the satellites in the geostationary orbit.It is more problematic for the satellites in the polar orbit because their position changes rapidly (some kilometres per second).The positional accuracy of an image has, however, always a major influence on results.

Accuracy of image matching and its influence on ACTH
The remaining sources of error in ACTH estimation are the accuracy of image matching and the accuracy of the wind velocity.Both have already been discussed by Seiz et al. (2007).
Their estimate of the height accuracy is based on the ratio between the base (a base is in photogrammetry defined as the distance between satellites projected to the Earth) and the height of the satellites.We did not follow their estimate because the heights of satellites used in the study are not the same (705 km for Aqua and Terra satellites carrying MODIS, and almost 36 000 km for MSG2 satellite carrying SEVIRI).
In addition, they consider the wind velocity that is in our approach already integrated by using two sequential SEVIRI images.Therefore, image matching is the only source of inaccuracy of the proposed method.It is more reliable if the compared images are made under similar conditions -for instance, the same acquisition time and the same geometry between the instruments and observing surface.If the observing surface changes while both images were acquired, image matching is likely to fail.This limits monitoring of the ash plume at the vent using the presented approach; 15 min is too long a period to observe the rising ash cloud at the vent because it may change its form and its height during this time interval significantly.Therefore, one might wonder why MSG1 data were not used in this study.MSG1 carries also SEVIRI just like MSG2, but since 2009 it operates in the rapid scan mode -it monitors only Europe at 5 min intervals.With this shorter retrieval interval, results would obviously be more robust.However, MSG1 is positioned approximately at 9 • E longitude and MSG2 at 0 • E. This 9 • difference does not seem much, but in the case of MSG2 Iceland appears already at the edge of the visible hemisphere.In the case of MSG1 data (located at 9 • E), Iceland is just on the horizon increasing the data spatial resolution and henceforth lowering ACTH accuracy.
In general the results of image matching between datasets acquired by two different instruments depend also on the used spectral bands.If the ash and meteorological clouds do not have similar reflectivity in MODIS and SEVIRI spectral bands, image matching is likely to fail too.Here we use the SEVIRI HRV band that covers a large part of visible and near infrared spectrum.Furthermore, we used MODIS band 1 that covers only the "red" spectrum of the visible light.The comparison of both datasets has shown that the data covering clouds and volcanic ash are highly correlated, justifying our bands' selection.In addition, different atmospheric path length of SEVIRI and MODIS measurements do not influence the results.For the sake of completeness, we note that image matching is able to detect also transparent clouds.This is usually not a problem if the underlying surface is homogenous.For instance, dark ocean surface is a suitable background.If the texture in the background varies significantly, the image matching will fail.
To analytically estimate the accuracy of image matching, we apply some simplifications: we ignore Earth's curvature (with the distance of the parallax being below 40 km this leads to a height error of not more than 100 m), assume that MODIS looks into nadir, and consider merely a 2-D space (meridian plane).These assumptions result in a simplified geometrical problem (Fig. 10 a); the error increases linearly with increasing error in parallax P ( ACTH = P • tan α).The satellite elevation angle α depends on the latitude (because of the simplification it actually depends on the angular distance from SEVIRI's nadir).Figure 10b shows the relationship between the ACTH error coefficient and the latitude.The ACTH error coefficient correlates the ACTH error with the coordinate error of the matching pixel being false by 1 km. Figure 10 shows that the error decreases with increasing latitude.In our case (assuming latitude of 60 • ) this coefficient equals to 0.4 km km −1 .If image matching results in a mistake of half of pixel to the North (this corresponds to 1.5 km in the north-south direction for the case study area), the ACTH error will be about 0.6 km.
The reliability of image matching can be estimated from the intersection distance between the MODIS and SEVIRI lines of sight (Eq.2).The reliability of the least-square solution of the lines intersection is large if the intersection distance between the lines of sight is close to zero.Mainly the north-south component is of interest in this case study area because the parallax between MODIS and SEVIRI extends approximately in this direction.The east-west component of the intersection distance has almost no influence on ACTH, thus the intersection distance was projected onto the north direction.These values are usually in the range of up to 1 km, which is less than half the pixel size for the study area.Furthermore, the intersection distance depends not only on the image matching accuracy, but also on the accuracy of the satellite position and positional accuracy of the image data.This makes the intersection distance a very suitable estimate for the general accuracy of the method.

Comparison with other methods
In order to further assess the validity of our approach we compare our results with ACTH presented in other studies.First, we compare our results with the plume height measured by MISR (Nelson et al., 2008).Wind-corrected heights of the ash plume were obtained from the MISR Plume Height Project (http://misr.jpl.nasa.gov/getData/accessData/MisrMinxPlumes/).Figure 11 shows the difference between our ACTH and MISR ACTH for 19 April at 12:51 UTC, which are significant.One reason for the discrepancy may be the vicinity to the vent; the northern part of the plume is probably still developing, and it moves also vertically, making the "measured" ACTH of MISR as well as ours quite unreliable.In the southern part of the plume the difference seem to be more systematic.Our "measurements" are about 1 km higher than the MISR ACTH.The reason for that most likely lies in the poor positional accuracy of the SEVIRI data.The visual comparison with MODIS data shows that Iceland (in SEVIRI data) is shifted by a pixel to the north.Such a shift was not observed in other parts of Europe (United Kingdom, or France).We conclude that the positional accuracy of SEVIRI data is not optimal with this inaccuracy propagating into ACTH (at least in the closest vicinity of Iceland).
We compare our results also with the standard MODIS cloud product MOD06 L2 based on the CO 2 absorption technique (Platnick et al., 2003) with a spatial resolution of about 5 km.The product estimates the cloud-top height pressure and not geometrical height.To convert the pressure into the geometrical height, we use atmospheric sounding data for the location of Tórshavn (Faroe Islands; 62.01 • N, 6.76 • W) on 15 April 2010 at 12:00 UTC.The time of the atmospheric sounding is not exactly the same as the time of MODIS retrieval (11:35 UTC; see Fig. 8a).The valid plume pixels seem in general to be higher in our results than in the MOD06 L2 product (Fig. 12a); the offset is approximately 650 m.
The same sounding data are used also to estimate ACTH with the BT method (Fig. 12b).Because the measured temperature in the sounding data fluctuates, it is possible to apply this sounding only for the clouds that are between 4 km and 11 km.The ash plume was certainly in this height range in the study area, but meteorological clouds around it most likely not.Thus Fig. 12 compares ACTH only in the area of the ash plume.The results of the BT method have an offset of about 2200 m compared to the values determined by our newly proposed method.An offset larger than 1000 m was already observed in similar previous studies (Genkova et al., 2007).The reason for this offset is related to the retrieval of a radiative height that does not estimate exactly the top of the cloud, but the height of the "cloud centre".

Conclusions
In our case study of the 2010 Eyjafjallajökull eruption, we applied the proposed method only to data collected in the visible spectrum, because of their high spatial resolution, but in general our technique it can also be used with e.g.infrared spectra.New instruments with better characteristics will be available in the years to come, allowing more elaborated applications of our method.For instance, the third generation Meteosat should provide also regional rapid scans every 2.5 min at a spatial resolution of 0.5 km in a visible band (0.6 µm) or 1.0 km in a thermal band (10.5 µm).Use of thermal bands will allow the application of the proposed method also during night-time.
In the future we plan to improve the accuracy of the image matching to the sub-pixel level.The accuracy of the current ACTH estimates is better than 0.6 km.With sub-pixel image matching, the estimates of shifts could be enhanced to about 10-20% of the pixel size.This would increase the ACTH accuracy to 0.2 km or even better.To achieve appropriate results in image matching data, pre-processing will require more effort: (1) the use of Wallis' filter to locally enhance contrast in the data, and (2) additional corrections of image registration by using vector data to improve the positional accuracy.
The greatest advantage of the proposed method to determine ACTH is its independence of physical assumptions.The method has a purely geometrical background, thus there is no need for any additional data (as for instance ash emissivity or atmospheric temperature profile).The geometrical solution is fast and it can be applied to any combination of a polar orbiting satellite and a geostationary satellite.The only requirement is that the temporal resolution of the geostationary satellite is high enough that clouds in the size of the used moving window do not change their form significantly, but just move to a new position.This allows correlating data retrieved from both orbits even if they were not retrieved at exactly the same time.The combination of geostationary and any polar orbiting instruments that have similar spectral bands can be used to measure ACTH a few times per day.It can be used as well for monitoring of meteorological clouds, volcanic plumes, mineral dust and vegetation fire plumes.

Fig. 1 .
Fig.1.If an ash cloud is observed simultaneously from two satellites the position of the cloud seems to be shifted for a parallax when comparing both datasets.The parallax is a consequence of viewing geometries of both satellites and the height of the ash cloud.

Fig. 2 .
Fig. 2. Aggregation of MODIS data into SEVIRI's grid system; the ash cloud is marked by the red oval.Top left: an example of original MODIS data (part of MODIS band 1 data showing the Eyjafjalljökull ash cloud over the North Sea on 15 April 2010 at 11:35 UTC).Top right: A scheme of SEVIRI's point spread function -each dot represents a MODIS measurement, thus it is coloured according to its position within the SEVIRI pixel (dark purple dots represent greater weight than pale yellow dots).Bottom left: MODIS data transformed in the SEVIRI grid system (see the difference in the overall geometry that is a consequence of a difference in viewing angles of both instruments).Bottom right: SEVIRI data retrieved at approximately the same time as MODIS (11:27 UTC).

Fig. 3 .
Fig. 3. Schema of area-based image matching.For more details see text.

Fig. 4 .
Fig. 4. Refinement of shifts using image pyramids.For more details see text.

Fig. 6 .
Fig. 6.BTD (data for 15 April at 20:00 UTC) reveals plumes dispersing at two different heights.Highly positive BTD values in dark blue over North Sea indicate ice crystals.Below it spreads a reddish plume with negative BTD.

Fig. 7 .
Fig. 7. SEVIRI data on 17 April at 13:12 UTC: (a) the cross-correlation index (see Eq. 1) between the SEVIRI image and MODIS image retrieved at 13:15 UTC; (b) optimal matching shift between both datasets in column direction; (c) optimal matching shift between both datasets in line direction; (d) parallax between both datasets; (e) intersection distance between lines of sight projected to the north direction; (f) ACTH for all available pixels; (g) ACTH for pixels where intersection distance projected to the north direction is lower than 1.5 km.Please note that b-d shows only CI and the shifts of the SEVIRI image retrieved before and not after the MODIS retrieval.

Fig. 10 .
Fig. 10.ACTH error coefficient as a function of the latitude.(a) shows the simplified geometry discussed in the text and (b) the variation of the error with latitude.

Fig. 12 .
Fig. 12. Differences between the valid ACTH of the ash plume and the standard MODIS top height product (a) and results of the BT method (b) on 15 April at 11:35 UTC (our ACTH are shown in Fig. 8a).