Ash fall impact on vegetation: a remote sensing approach of the Oldoinyo Lengai 2007–08 eruption

Abstract

Impacts from ash fallout on the environment can be widespread and long lasting, even from moderate-size eruptions. Assessing ash impact on vegetation and the indirect impacts for people is often difficult in the field. Here it is assessed how satellite data can help to map vegetation affected by ash and how temporal analysis enables characterization of vegetation recovery rate. The 2007–08 eruption of Oldoinyo Lengai, north Tanzania, is here used as a case study. An 8 year-long (2005–2012) time series of half-monthly average of the Normalized Differential Vegetation Index (NDVI) is constructed at 250 m spatial resolution from the Moderate Resolution Image Spectro-radiometer (MODIS) sensor. Interpolated rainfall data is used to isolate NDVI values departing from the normal seasonal cycles. Month-to-month NDVI comparison, linear temporal trend analysis and Principal Component Analysis enable to identify a 11 × 4 km area over which ash fallout significantly affected the state of the vegetation. After the eruption’s end, time series of various recovery indices highlight a circumferential pattern in vegetation recovery. The estimated recovery time varies from more than 5 years to less than 6 months with increasing distance from the volcano. A non-linear moderate, but statistically significant, relationship is found between the recovery indices and the spatial variation of ash thicknesses measured in the field. Combining field and remote sensing constraints enable to re-assess the volume of the eruption to ~2 × 10⁷ m³. The spatial pattern of the ash-affected area matches with the spatial contrast in the impact experienced by the local communities. The method applied here opens the scope to document impact and intensity of ash fallout in areas where systematic field work is not possible and to support recovery plans for populations affected by ash fallout.

Introduction

Ash fallout produced by volcanic eruptions, even of moderate magnitude, are responsible for severe impacts around volcanoes, due to burial, overloading, or different physical or chemical processes affecting plant growth or foliage state even for thin tephra coverage (see Ayiris and Delmelle 2012 and reference therein). Fires ignited by hot tephra and chemical impacts of toxic gases via air pollution, acid rains or soil water contamination are other processes by which explosive eruptions impact on vegetation (e.g., Delmelle 2003; Delmelle et al. 2002). Tephra layer thickness, grainsize and physico-chemical properties as well as the morphology and phenology of the plants and the residence time of the tephra layer have been shown to impact the extent of these impacts (Ayiris and Delmelle 2012). Vegetation destruction or deterioration by volcanic ash may result in long term impacts for livestock, agriculture and eventually for human activity and health (e.g., Wilson et al. 2011). Shrub and herbaceous vegetation have been shown to be more vulnerable than trees to ash, being already severely affected by a few centimetres of ash (del Moral and Grishin 1999; Ayiris and Delmelle 2012).

It is generally time-consuming and difficult to accurately constrain the extent and thickness of volcanic ash layers, and the spatial variation of impacts on the flora, fauna and human activities from field surveys. This is especially true for long lasting eruptions when ash fallout will occur repetitively over weeks or months, for volcanic regions which are less accessible and for more distal and thin ash layers that are poorly preserved over more than a season.

Although many studies highlight the impact of volcanic ash fallout on vegetation (e.g., Ayiris and Delmelle 2012 and references therein), fewer studies document the recovery of vegetation after ash deposition and the factors influencing its recovery, with the exception of the studies following the eruptions of Mt Lamington (Taylor 1957), Pinatubo (De Rose et al. 2011), Mt St Helens (Saito et al. 1998; Frenzen 1992; Harrington et al. 1998; Marzen et al. 2011,) and Chaiten volcanoes (Martin et al. 2009). Vegetation types, ash thickness, rainfall, erosion pattern and rate, and human influence were identified as important factors controlling vegetation recovery (Frenzen 1992; Harrington et al. 1998; del Moral and Grishin 1999; Marzen et al. 2011). De Rose et al. (2011) further concluded that watershed topography and geology have an important impact on the trajectory and the speed of vegetation recovery.

Most of these studies were based on temporal analysis of species distribution documented in the field. Only De Rose et al. (2011) and Marzen et al. (2011) used vegetation indices derived from remote sensing data to analyse vegetation recovery at volcanoes. Marzen et al. (2011) monitored vegetation recovery from 1980 to 2005 at Mt St Helens using bi-temporal difference analysis of NDVI-Landsat images. Such an approach is more commonly used to map the impact of, and study the recovery from, forest fires (Diaz-Delgado et al. 1998; Salvador and Pons 1996; Riaño et al. 2002; Diaz-Delgado et al. 2003; Gouveia et al. 2010; Hope et al. 2012).

Here we study the temporal evolution of vegetation around Oldoinyo Lengai (OL) in North Tanzania using low spatial resolution remote sensing data. The objective is to investigate how the vegetation characteristics responded to the 10 month-long violent Strombolian to Vulcanian-style eruption in 2007–08 (Vaughan et al. 2008; Kervyn et al. 2010). It is tested if the analysis of vegetation change from remote sensing time series can help to identify and map the zone that has been affected by significant ash fallout. In addition, using the analogy with vegetation recovery after forest fires (Diaz-Delgado et al. 1998; Salvador and Pons 1996; Riaño et al. 2002; Diaz-Delgado et al. 2003; Gouveia et al. 2010; Hope et al. 2012), we investigate how the vegetation recovery rate can be quantified and if it is correlated with the spatial variation of the thickness of the deposited ash layer.

Oldoinyo Lengai and the 2007–08 eruption

Oldoinyo Lengai (OL) is a stratovolcano of the East African Rift System, located in North Tanzania, 16 km south of Lake Natron (Fig. 1a). With its 2952 m a.s.l summit, Oldoinyo Lengai rises 2000 m above the plain of the Gregory rift valley. Similar to the Serengeti plains to the West, the vegetation in the rift and above the rift escarpment is characterized by open grasslands to shrublands with scattered trees, varying with topography and rainfall distribution (Fig. 2f; see Reed et al. 2009; White 1986). The Crater Highlands south west of OL are covered by a dense shrubbed forest (Fig. 1b; Reed et al. 2009). The climate of the region is characterized by a short rain season from October to December, and a long and more abundant one from March to May.

Fig. 1

Oldoinyo Lengai, surrounding volcanoes and tectonic structures. a: Shaded relief map based on the SRTM Digital Elevation Model highlighting main volcanoes and craters; (b) Land cover map of the study area based on the FAO classification scheme (FAO 1998). Red lines highlight the approximated isopach lines of the 2007–08 ash layers inferred by Sherrod et al. (2013). The full triangle marks the location of Oldoinyo Lengai

Fig. 2

Illustration of the impact of ash on Oldoinyo Lengai and its direct surrounding. a: Summit of OL looking south-east during the eruption (March 12, 2008; photo courtesy B. Wihelmi); (b) Ash covering vegetation along the west flank of OL and the rift escarpment (January 18, 2008; photo courtesy T. Pfeiffer); (c) Ash fall in valley cutting across the rift escarpment two years after the eruption (June 8, 2010; photo courtesy Th. Holden); (d) Boma Lesele, 4 km west of OL crater located in the major ashfall zone (December 6, 2007; photo courtesy of J. Keller); (e) Ash fall deposits on W flank of OL and destroyed vegetation (May 14, 2011; photo by S. Stadlin); (f) Recent regeneration of vegetation north of OL (May 2013; photo by S. Stadlin)

Oldoinyo Lengai has been extensively studied for its petrology, as it is the only active volcano known to have erupted natrocarbonatite lava from 1983 to 2006 (e.g., Bell et al. 1998; Dawson 1998; Zaitsev and Keller 2006; Kervyn et al. 2008). OL is however built up dominantly by nephelinite and phonolite ash and lava deposits, natrocarbonatites representing a minor volume (e.g., Klaudius and Keller 2006). OL eruptive activity is characterized by phases of effusive eruption of natrocarbonatite lava and nephelinite explosive phases, known to have occurred in 1917, 1940–1941, and 1966–1967 (Dawson et al. 1968; 1995; 1992).

On September 4, 2007, a 3 km-high eruption column dispersed ash up to 20 km from the summit (Vaughan et al. 2008). This event marked the start of a new explosive eruptive phase characterized by repetitive short-lived ash eruptions of varying intensity, as documented by Kervyn et al. (2010). After the onset, the activity varied from 100 m high ash jets to 2–15 km high, steady or unsteady, eruption columns. The eruption was most intense at the end of September 2007 and in February-March 2008, the latter phase being associated with the largest ash dispersal and the first pyroclastic flows to be documented at OL (Kervyn et al. 2010). The petrology of the 2007–08 eruption and the processes responsible for the eruption style transition are discussed in Mitchell and Dawson (2007); Keller et al. (2010); de Moor et al. (2013) and Bosshard-Stadlin et al. (2014). The main explosive activity lasted until end of April 2008 with only gas and minor ash emissions being observed afterward until resumption of lava effusion at the end of 2008. The volume erupted by the 2007–08 eruption sequence was estimated at ~10⁷ m³ (Kervyn et al. 2010) but this estimate was not based on systematic documentation of ash thickness in the field.

Most of the erupted material was deposited on the upper volcano flanks. Proximal deposition built up a ~400 m wide, ~75 m high pyroclastic cone. By early December 2007, Keller et al. (2010) observed 5 and 1 cm ash thickness at the W and S volcano base, respectively. Three months after the onset of the eruption, all the vegetation on OL’s flanks was already dead or fully buried (Fig. 2). Even if ash fallout was reported in villages within the rift valley for specific eruptive events (e.g., in Engare Sero, north of OL, on September 4, 2007), most of the fallout was recorded on the W flank and above the 500 m high rift escarpment. This is consistent with the West-blowing trade winds (Kervyn et al. 2010). The most intense eruptions in February-March 2008 dispersed ash over larger areas, with millimetre-thick fallout being reported as far as 100 km to the SW (i.e., Lake Eyasi and Ngorongoro caldera; Kervyn et al. 2010).

A preliminary isopach map of the 2007–08 ash fallout was produced by Sherrod et al. (Sherrod et al. 2013), constraining a 5–10 cm isopach line, and partially constraining the limit of a 1–5 cm thickness isopach line (Fig. 1b). Naiyobi and Kapenjiro villages, each with 5000 to 6000 inhabitants, were affected by the thickest ashfall of all settlements surrounding OL. This first characterization of the ash fallout is here completed with additional field observations and remote sensing analysis.

No study exists on the processes by which OL eruption impacted the surrounding vegetation but field observations (Fig. 2) suggest that burial, overloading or prolonged coverage of leaves were the main processes causing vegetation destruction or decay. Fires ignited by lava overflowing the crater in the days before the onset of explosive eruptions also impacted vegetation on OL west and east flanks (Kervyn et al. 2010). Gas emissions from OL have been shown to be dominated by water vapour and CO₂, with limited emission of SO₂, H₂S, HCl or HF (Koepenick et al. 1996; de Moor et al. 2013). Direct effect of gas emissions from OL has never been reported and is considered negligible here. Ash leaching experiments carried out by de Moor et al. (2008) and Stadlin (2015) showed that OL ash contained a high concentration of leachable ions, mostly fluorine, which constitutes a major hazard for local groundwater supply and for the grazing cattle but without observed direct effect on vegetation state.

The local communities living around OL are dominantly Masai who are dependent on cattle for their livelihood, but other tribes developing agriculture are also present in the region. Ash fallout was reported to cause respiratory problems to local residents living in villages on the escarpment west of OL. At locations where heavy ash fallout was experienced, within 25 km west of the volcano, local communities spontaneously evacuated to less affected areas. A survey conducted in May 2011 in 8 different settlements (99 interviewees) surrounding OL highlights the difference in impacts experienced based on the proximity and geographic location of the settlement to OL. In Naiyobi, a large settlement located downwind from OL, over 80 % of interviewees decided to relocate on a voluntary basis for several months during the eruption. Temporal relocation of up to 10 months also occurred for part of the population in Engare Sero and Engaruka-Ngwesuku, 15 and 7 km of OL in the N and SE direction, respectively. In the other settlements located further away from OL, only few respondents declared to have evacuated, typically for shorter durations.

The decision of relocation was motivated mainly by the lack of green pastures for the cattle (NEMC 2008), but interviewees were also worried about their safety and health and access to safe water. Most of them reported eye irritation and some reported respiratory problems. The impact for the cattle included death, miscarriage and skin irritation. 74 % of respondents experienced crop destruction due to the ash and 88 % reported an overall loss of income. Thin ash fallout was enough to prevent livestock from grazing and forced temporary resettlements. The total number of displaced people is evaluated at 5,000 to 7,000 people, whereas about 65,000 people are estimated to have been affected by the eruption (NEMC 2008). Several NGO’s and the Red Cross had to provide food aid for approximately 36,000 people (Msami 2007).

Datasets and methods

Vegetation index: dataset and pre-processing

Different Vegetation Indices (VI) have been defined to characterize the state and density of vegetation from the spectral reflectance of a surface in the visible and near-infrared wavelengths. The Normalized Difference Vegetation Index (NDVI), defined by the difference of reflectance in near-infrared and red, divided by the sum of both, was one of the first indices proposed (Tucker 1979) and up till today the most used. Different Soil Adjusted VIs (SAVI) were defined and tested to better relate the VI to the biophysical properties of the vegetation but these indices do not always outperform the NDVI (Carreiras et al. 2006; Clemente et al. 2009). Veraverbeke et al. (2012) concluded that in applications concerning post-fire vegetation recovery, the NDVI is the most suitable index due to its normalizing capacity for vegetation variability, which makes that it is strongly correlated to the above-ground biomass in different ecosystems (Clemente et al. 2009; Lhermitte et al. 2010). NDVI will also be used here as it has proved robust in previous remote sensing studies on vegetation monitoring, including vegetation recovery after forest fire (e.g., Diaz-Delgado et al. 1998; Riaño et al. 2002; Hope et al. 2012).

To monitor vegetation evolution during and after ash fallout, a dataset offering a balance between a sufficiently high spatial and a moderate temporal resolution was selected. In addition the dataset had to be able to provide reliable and consistent data throughout the year despite the frequent cloud coverage common to most equatorial countries. The Moderate Resolution Imaging Spectro-radiometer (MODIS) satellite acquires daily images of the entire world at 250–1000 m spatial resolution since 1999 (Ernst et al. 2008). The MOD13Q1 16-day NDVI composites are used here (LPDAAC 2013; http://lpdaac.usgs.gov): these data are derived from the MOD09 daily surface reflectance product, which is already corrected for molecular scattering, ozone absorption and aerosols (Vermote et al. 2002). These composite images are built up using the Maximum Value Composite algorithm filtering for data quality, cloud and viewing geometry (Huete et al. 2002; Solano et al. 2010). The first 16-day composites of each month from January 2005 till May 2013 were used to build up a time series for the study area.

As errors due to cloud cover, sensor characteristics and variations in atmospheric conditions are still present in the MODIS NDVI time series (Viovy et al. 1992), each pixel’s time series was filtered with the TIMESAT algorithm. TIMESAT (Jönsson and Eklundh 2002, 2004) aims at least squares fitting of a function on the time series to filter out noise and remove spikes (Olsson et al. 2012). Contrasted fitting functions, Gaussian or logistic, were obtained for different land-cover categories. The result of this operation was the corrected monthly NDVI_corr time series for each pixel over 7 years (Fig. 3).

Fig. 3

Typical NDVI time series for selected pixels in the study area: (a) on Oldoinyo Lengai, (b) in the ash affected zone, (c) in the Crater Highland forest, (d) in the savannah. The original NDVI and the NDVI_corr after TIMESAT correction are shown for each pixel from May 2005 to May 2013. (e) Interpolated rainfall time series from May 2005 to May 2010. Note the grey shading highlighting the eruption period. f: NDVI image of study area in February 2008 with location of pixels analysed in A-D

Seasonal variation in precipitation is the main control on the temporal variation of the NDVI_corr time series (Ding et al. 2007; Markogianni et al. 2012). In order to isolate the effect of ash on vegetation, the NDVI time series had to be de-correlated from precipitation (Zhou et al. 2009). Daily precipitation data were acquired from the Tanzania Meteorological Agency from March 2005 to May 2010 for 4 stations located within 110 km of Oldoinyo Lengai (Arusha 110 km to SE, Monduli 88 km to SE, Loliondo 80 km to NNW, Karatu 70 km to SW). All stations showed similar seasonal and inter-annual rainfall variations, despite different total annual precipitation. A single rainfall time series was interpolated based on an Inverse Squared Distance Weighted interpolation from the 4 stations. The precipitation value obtained for the location of Oldoinyo Lengai was used to decorrelate the NDVI_corr time series over the entire study area. This simplification is required due to the scarcity of available weather stations and the lack of local modelling of the rainfall spatial variability in the study area.

Zhou et al. (2009) used regression analysis to filter out NDVI variation caused by rainfall. Here we divided the study area into sub-regions of different land cover and elevation ranges, assuming that the reaction of vegetation to precipitation depends on altitude and vegetation type (Reed et al. 2009). Rainfall decorrelation was carried out separately for each sub-region using a polynomial relationship between NDVI and preceding precipitation, similarly to Razali and Nuruddin (2011). As the ash-affected vegetation would disturb the NDVI-precipitation relationship, the parameters of the regression for each sub-region were obtained using only pixels located outside a broadly defined zone potentially affected by ash.

Precipitation amounts were summed up for every 16-day period to correspond to the time span covered by the semi-monthly NDVI data. The second-order regression analysis was carried out between the precipitation of the preceding 16 day-period and the corresponding average NDVI value for each vegetation-elevation class separately. This procedure was repeated for preceding rainfall periods of 32 to 102 days to obtain the time lag which showed the highest correlation between NDVI and precipitation. NDVI variation expected from rainfall variability (NDVI _exp ) can be derived with the following function:

$$ NDV{I}_{exp} = a + b\ .\ {P}_x + c\ .\ {P_x}^2 $$

(1)

with a, b and c the coefficients of the regression analysis and P_x the preceding precipitation for the period returning the highest correlation (Fig. 4, Zhou et al. 2009). This way, the deviance between the expected and observed value of NDVI_corr can be derived:

Fig. 4

Example of the polynomial relationship between the average bi-monthly NDVI and the 48 days antecedent rainfall for herbaceous vegetation between 1500 and 2100 m a.s.l. The 48 days antecedent rainfall provided the best statistical relationship

$$ NDV{I}_{dif} = NDV{I}_{corr}\hbox{--}\ NDV{I}_{exp} $$

(2)

A negative NDVI_dif means a lower NDVI is observed for a pixel than the value expected for that vegetation type at that elevation after the experienced rainfall. It will here be tested if the variation of NDVI_dif can be attributed to the impact of ash fallout.

Mapping vegetation affected by ash

Detection of the ash impact on vegetation around OL was performed using bi-temporal analysis and Principal Component Analysis of the NDVI_dif time series. Bi-temporal change detection compares two periods in time by removing all variances caused by non-relevant factors to measure change only caused by differences in the variable of interest (Green et al. 1994). Using bi-temporal change detection, the sensor’s spectral and spatial resolution, the acquisition date and time, the atmospheric conditions, and the phenological state of the vegetation should be similar (Lu et al. 2004). The bi-temporal change analysis was carried out on the NDVI_dif time series, removing the potential effect of rainfall variation. Data were compared for the same month of the year for the periods preceding, spanning and following the eruption, after normalization of the data range in each scene.

A Principal Component Analysis (PCA) linearly transforms data into a new set of uncorrelated principal components sorted according to the variance explained in the original data (Schowengerdt 2007). The standardized PCA algorithm converts all scenes of NDVI_dif values to standard scores which results in equal weights for all scenes within the time series in the derivation of new component images (Schowengerdt 2007). A standardized PCA was carried out on the NDVI_dif time series from July 2005 to May 2010 with the objective of isolating the signal of ash-affected vegetation. The duration of the time series was chosen in order to cover two full years before and after the eruption.

Vegetation recovery indices

In addition to mapping the area affected by ash fall with remote sensing, we aim at characterizing the recovery rate of the affected vegetation and see how this recovery is related to the ash thickness deposited at any given location. In order to be able to examine a long enough time period for deriving Vegetation Recovery Indexes, the NDVI_corr time series was extended to May 2013. NDVI_corr is used here, and not NDVI_dif, as the goal is to see how vegetation recovers to its pre-eruption state, including its sensitivity to rainfall variation.

First, an NDVI trend analysis was carried out using the Ordinary Least Square procedure (OLS). OLS calculates the slope coefficient α of the linear regression between time and values of NDVI_corr.

$$ NDV{I}_{corr}(t) = a\ .\ t + b $$

(3)

This results in an expression for the rate of change over a set of time steps in such a way that the slope values calculated for the period following the eruption can be interpreted as recovery rates. In a normal phenological year the beginning of a wet season is characterized by positive slope values, followed by negative slope values in the dry season leading to a slope value around zero for the yearly variation. This makes it possible to detect a disturbance which will cause these slope values to differ from zero. The NDVI trend analysis was carried out for three time spans: the pre-eruptive period (May 2005 – May 2007), the eruption (May 2007 – May 2008), the post-eruptive period as a whole (May 2008 – May 2013).

Second, different recovery indices were calculated for the different time intervals following the eruption (Table 1). Previous studies proposed vegetation recovery indices based on spatial comparison between undisturbed and deteriorated vegetation within the same scene (e.g., Diaz-Delgado et al. 1998; Diaz-Delgado et al. 2003; Riaño et al. 2002) or through temporal comparison of the vegetation state within each pixel separately (e.g., Gouveia et al. 2010; Hope et al. 2012; Lu et al. 2012). The high vegetation variability within the study area makes it difficult to select unambiguous control sites for spatial comparison. We therefore opted for the latter approach. As highlighted in Table 1, all the indices that were used are based on the pixel NDVI_corr value at different times t after the eruption (NDVI _post ) and the corresponding NDVI_corr value before (NDVI _pre ) and/or directly after the eruption (NDVI _erup or NDVI _lowest ).

Table 1 Vegetation recovery indices used to document the recovery of vegetation over time after the eruption. All NDVI values used in these indices are NDVI_corr values

The ‘Lack of greenness’ index (LG) is a simple difference between the NDVI value of each pixel before and after the event (Gouveia et al. 2010). The pre-event value is determined as the average NDVI for the corresponding month recorded from 2005 to September 2007 (NDVI* _pre ). The ‘Change in NDVI’ parameter (T_c) is similar to the ‘Lack of greenness’, with opposite sign, but it is normalized by the value of the pre-event NDVI (Lu et al. 2012), whereas the ‘Stand Regeneration Index’ (SRI) is a simple ratio between the post- and pre-event NDVI. The ‘Vegetation Recovery Rate’ (VRR) is a ratio of the recovery after the eruption and the change in vegetation caused by the eruption: the NDVI value recorded for the month following the end of the eruption (NDVI_erup) is subtracted from the pre- and post-event NDVI values (Chou et al. 2009; Lu et al. 2012). Finally, we propose the Ash Recovery Index (ARI) which is an adaptation of the VRR index based on Hope et al. (2012). In the ARI, the pre-eruption NDVI is calculated as the maximum NDVI value recorded for each pixel in the pre-eruptive period. The lowest NDVI (NDVI _lowest ) recorded in the year following the end of the eruption is used as representative of the vegetation status following the eruption.

In order to remove signals from pixels not affected by the eruption, the above indices were only calculated for the pixels for which the NDVI_{max pre} – NDVI _erup value is above 0.25. It should be noted that the value of LG will be negative and that of T_c positive but both decrease in absolute value over time for vegetation recovering after the eruption. The three other indices (SRI, VRR ARI) are positive and increasing over time for recovering vegetation.

Vegetation recovery rate

In order to define the recovery rate from the recovery indices (RI), a logarithmic regression model is fitted through the values of the RI through time, similar to Diaz-Delgado et al. (1998 1998).

$$ RI\ (t) = \alpha + \beta log\ t $$

(4)

With α the constant and β the slope of the logarithmic relationship, and t the time of the NDVI_corr image. This relationship is obtained for each pixel of the study area separately. The logarithmic relationship assumes that the recovery rate will be rapid directly after the event and will thereafter slow down to reach the pre-event level. It is expected that the slope of this relationship will give an indication about the intensity of the vegetation recovery and can thus be correlated to ash fall thickness. The type and strength of the relationship with ash thickness will be tested statistically. Negative slope values indicate no significant recovery through time and are therefore disregarded in the analysis. However, in order for pixels which have not yet recovered to be taken into account, a small constant is added to the derived slope values, in order for slope values close to 0 to become positive. This is justified as only pixels very close to the volcano display near-zero values, whereas vegetation not significantly affected by the eruption displays significant negative slope values (see Results).

Using this logarithmic regression it is also possible to define the recovery time for each pixel. The average value of the recovery index is calculated for NDVI_corr values observed during the pre-eruptive phase (January 2005 – June 2007). The recovery time is considered as the time required for the post-eruption recovery index to reach the mean value of the index in the pre-eruptive phase.

Validation dataset: field mapping of ash thickness

Thickness of the ash layer deposited by the 2007–08 eruption of Oldoinyo Lengai was documented in the field by several research parties. 20 observations were made in 2009 and 2010 by D. Sherrod (USGS, pers. comm.). de Moor et al. (2013) took 19 control points with ash thicknesses as averages of ~5-6 holes per site in 2009. An additional 105 locations with ash thickness observations, presented in Bosshard-Stadlin et al. (2014), were collected during field work between 2010 and 2013. It is to be noted that most sites are located within 7 km of OL summit, and that no observations have been made above the rift escarpment west of the volcano due to limited accessibility. This prevents the characterization of the medial to distal portion of the ash fallout. Recorded ash thickness varies between 4 mm and 84 cm, all thicknesses greater than 20 cm being recorded on OL flanks.

Results

Mapping ash impact on vegetation

Bi-temporal changes

Figure 5 illustrates the results from bi-temporal change detection. Limited vegetation changes are observed between October 2006 and October 2007 around OL, one month after the start of the eruption, beyond the general decrease in vegetation status over the area at the end of the dry season (Fig. 5a). Larger NDVI differences are obtained for the January 2007–2008 period over the entire study area (Fig. 5b). A 7 × 18 km area extending WNW from OL displays a significant negative vegetation change. January is in the middle of the growing season leading to greater differences in NDVI_dif also outside the ash-affected areas: the forest and shrub savannah region above the rift escarpment show a general improvement whereas the grassland savannah of the rift valley show a general decrease over that period, probably caused by local variations of the rainfall amount. Comparing the vegetation before and at the end of the eruption in May 2007 and 2008, the vegetation affected by the ash deposits is clearly identified above the background variations (Fig. 5c). The ash affected area is mapped as a near-elliptical 11 × 24 km area elongated to the west.

Fig. 5

Bi-temporal change detection results for NDVI_dif for different periods: (a) Oct 2006–07; (b) Jan 2007–2008; (c) May 2007–2008

Principal component analysis

Table 2 and Fig. 6 illustrate the results of the Principal Component Analysis on the monthly NDVI_dif time series from July 2005 to May 2010. The first six principal components have eigenvalues above one and explain together more than 80 % of the total variance (Table 2). PC1 reflects the main variation of the data (~60 % of the total variance). PC 2 to PC 5 display various spatial patterns and temporal variation in the loadings which can be interpreted as spatial and temporal variations in vegetation state not properly accounted for by the NDVI rainfall correction.

Table 2 Statistical results of the Principal Component Analysis, with the number of principal components considered in the interpretation marked in bold

Fig. 6

Principal component analysis of the vegetation time series: (a) Map of the factor scores for the sixth component of the PCA carried out on the NDVI_dif time series from July 2005 to May 2010; (b) Component loadings for the sixth principal component, with variables being the monthly NDVI_dif values for each pixel. Note the grey shading highlighting the eruption period

PC 6 explains 1.79 % of the total dataset. This component clearly isolates the signal caused by the ash fallout on the vegetation surrounding OL. Whereas almost the entire study area displays slightly positive values, the ash affected zone is characterized by highly negative values. This interpretation is consistent with the positive loadings of the months at the end of the eruption which corresponds also to the peak of the 2008 rain season. It is at the end of the growth season (i.e., April-May) that the deviation of the NDVI relative to its expected values based on the precipitation is the highest in the ash fall-affected area.

Vegetation recovery

Trend analysis

Figure 7 illustrates the results for the NDVI linear trend analysis on NDVI_corr time series. The pre-eruption slope values are generally close to zero, with stable vegetation status for the forested areas and more positive values for shrub and herbaceous vegetation, due to a longer and more abundant rain season in 2007 compared to 2005 and 2006. A larger range in slope values is obtained during the eruption (Fig. 7b). Most of the study area displays a positive slope, due to the continuous vegetation improvement associated with an abundant rain season in early 2008. Forested areas display a low sensitivity to precipitation with limited inter-annual variation of NDVI. The zone surrounding Oldoinyo Lengai is characterized by highly negative slope values, due to the rapid degradation of the vegetation state impacted by ash deposition. The steepness of this decreasing trend varies with distance from the volcano summit and the dispersion axis (Fig. 7b).

Fig. 7

Slope values resulting from the OLS trend analysis. Analysis is carried out on NDVI_corr values per pixel for three different time periods: (a) pre-eruptive phase May 2005 – May 2007; (b) during the year of the eruption May 2007- May 2008; (c) post-eruptive phase May 2008 - May 2013; (d) Schematic representation of the temporal evolution of the NDVI for an area affected by ash and the associated trend lines for the three period of interest shown in a, b and c

In the post-eruptive phase from May 2008-May 2013, the highest positive slope values are observed in the outer part of the ash fall zone as defined by the principal component analysis (Fig. 7c). These pixels are recovering, and have a higher increase in NDVI during the growing season compared to areas that have not been affected by ash deposition (Fig. 7d). In the direct proximity of OL, the slope values are close to zero suggesting limited to no recovery of the vegetation. It should be noted that the range of slope values is smaller than in previous time periods, due to the longer duration of the time series.

Vegetation recovery index

The transformation of NDVI_corr values to recovery indices is carried out to highlight the evolution of the ash-affected vegetation after the eruption ends. The different recovery indexes compare, on a pixel-to-pixel basis, the vegetation state after the eruption with its state before the eruption (see section 3.3; Table 1).

Figure 8a illustrates the Lack of Greenness (LG), which simply subtracts the pre-eruption value from the post-eruptive value, 1, 12 and 36 months after the end of the eruption. Directly after the eruption, the ash-affected area is characterized by strongly negative LG values, with values increasing towards the edge of the affected zone. Up to 25 km West of OL, the LG values are significantly below the variation observed within the surrounding area, suggesting a significant impact of ash on vegetation. After 3 years, the distal parts of the ash affected zone have totally recovered, with a vegetation status better than before the eruption (positive LG values). Only a limited zone extending 10 and 7 km west of OL still displays a significant impact in May 2009 and 2011, respectively.

Fig. 8

Vegetation recovery indices for the ash affected area and its direct surroundings 1 month (first row), 12 months (second row) and 36 months (third row) after the end of the eruption. The graphs at the bottom give the temporal evolution of the recovery index value for the pixel marked by a red point in A for the period May 2008 to May 2013. a: Lack of Greenness (LG); (b) Change in NDVI (T_c); (c) Ash Recovery Index (ARI)

Figure 8b illustrates the change in NDVI (T_c). T_c returns a high value for ash-affected vegetation and its value decreases with vegetation recovery. A similar spatial pattern is observed as for the LG index but the contrast within the ash-affected area is higher.

Figure 8c shows the Ash Recovery Index (ARI). Low ARI values characterize the ash fall zone whereas vegetation outside the ash-affected zone shows values close to, or above, one. Although the recovery dynamics of the ash-affected area is well characterized by this index, the contrast is low with the vegetation not-affected by the ash, in the north-west and south-east of the study area. After 3 years, the values of the ARI index in the central part of the ash-affected zone are similar to the other parts of the study area, but the edge of the ash-affected zone displays higher values.

All three indices succeed in highlighting the spatio-temporal evolution of the vegetation recovery within the ash-affected area. The index values one month after the eruption match the mapping of the ash-affected vegetation with PCA, but extent further to the West (Figs. 5, 6 and 7). With time, the area affected by thinner ash deposition recovers, and the zone where vegetation remains affected diminishes in size. After 3 years, only a 3 x 7 km area still shows degraded vegetation. Also, the edge of the ash-affected area displays high values of the indices (mostly with LG and ARI): this suggests that vegetation which recovers after being affected by a thin ash cover is in a better state than before the eruption.

Logarithmic regression

The slope of the logarithmic regression of the recovery indices through time reflects the recovery rate. Half and full recovery time were calculated as the time required for each pixel to reach half or the same value of the recovery index respectively, as calculated before the eruption.

Figure 9 shows the slope of the logarithmic regression and the derived recovery time per pixel for the ARI index. All recovery indices resulted in the same spatial pattern for parameters of the logarithmic regression fitting. The ash-affected area surrounding OL can be sub-divided into two zones (Fig. 9a). The proximal area, limited to OL west flank and basis affected by the thickest ash deposition, is characterized by low to zero slope values. These indicate that there is a very slow recovery or no recovery. The external part of the ash-affected zone in the main wind direction is characterized by high slope values. These indicate that there is a rapid recovery. One to three years is estimated as the time for full recovery in these peripheral zones of the ash-fall area. For the central part, the recovery was not completed after 5 years (Fig. 9b and c).

Fig. 9

Map of the logarithmic regression parameters for fitting the Ash Recovery Index (ARI) values through time. (a) Slope values, (b) Half recovery time (c) Full recovery time. (d) Logarithmic regression between the measured ash thickness and the slope values of the ARI. (e) Spatial distribution of field points for which thickness is documented. Grey values in D and E marks points for which no positive vegetation evolution is observed. The purple ellipse in a, b, c and e show the location of the 3 cm isopach line used for volume interpolation

Relationship with ash thickness

Also examined was whether the recovery rate of the vegetation, measured by the slope of the logarithmic regression, can be related to the thickness of ash measured in the field. Using 144 ash thickness points measured in the field, all recovery indices showed a general power law relationship with the observed ash thickness. The statistically best fit is found for the Ash Recovery Index (R² = 0.28), followed by the Change in NDVI Tc (R² = 0.25) with weaker fitting for other indices. Points where a thin ash layer was observed are characterized by high value of the recovery slope whereas a thicker ash layer corresponds to a lower value of slope recovery indices. In addition to the limited spatial and thickness ranges of field data, the large scatter around these power law trends highlights the fact that ash thickness is only but one factor affecting the recovery rate and that the relationship between ash thickness and recovery rate might not be a continuous one but could be controlled by ash thickness thresholds. Differences in the statistical relationships between the different indices might be attributed to the capacity of the various indices in capturing the dynamics of the vegetation recovery and their sensitivity to seasonal and inter-annual rainfall variations.

Some field measured ash thicknesses do not fit the relationship as they are characterized by a negative slope for the logarithmic regression model, i.e., the vegetation state is not improving in the five years following the eruption (Fig. 9d). Negative recovery slopes are obtained for points characterized by a thin ash layer. Of the 15 to 22 points (10–15 % of total dataset) that had to be excluded from the analysis, two were considered as abnormal measurements: over 20 cm of ash was measured by Sherrod more than 7 km north-east of the volcano. These observations do not match with the other points and could be attributed to local ash reworking and accumulation or erroneous identification of the layer caused by the 2007–08 eruption. Over 80 % of the remaining excluded points (10–13 data points) have ash thicknesses of 3 cm or less and are located upwind from the volcano (Fig. 9e). 30 % of the points that do display a significant vegetation recovery are also characterized by ash thicknesses of 3 cm or less. However, it can be considered that 3 cm of ash is the threshold thickness above which all vegetation types in the study area showed a significant recovery trend after the eruption, and thus the ash thickness required to systematically affect the state of vegetation.

Eruption volume estimation

The extent of the ash-affected vegetation can be used to estimate the eruptive volume of the OL 2007–2008 eruption. As the field data are insufficient to define a detailed deposit thinning trend based on a sufficient number of isopach lines (e.g., Bonadonna and Costa 2012), we used a spatial interpolation technique to derive the thickness distribution based on available constraints. In addition to the field-measured ash-thickness points, the following assumptions were made to constrain the ash volume: 1. the largest field-measured ash thickness (84 cm) was also attributed to OL summit; 2. the elliptical zone of ash-affected vegetation identified in Fig. 9 was defined as the 3 cm isopach, based on the above defined threshold; 3. two additional ellipses reaching 50 and 100 km downwind from OL where defined for the 1 and 0 mm ash thickness, respectively. The latter constraints are based on the largest distance of reports mentioning ash fallout. These boundaries are poorly defined but are needed to constrain the volume: their assumed position is a conservative one. Based on a total of 705 points, including field observations and sampling along the defined isopach lines, the thickness of the ash layer was interpolated using Inverse Squared Distance Weighted Interpolation (Fig. 10). Although the exact shape of the ash fallout area is not represented by this interpolation, it does account for the observed thickness decrease with distance. Integrating this ash fall thickness map gave a new eruptive volume of 1.52 10⁷ m³. 43 % of this volume is deposited on OL flanks, within 3.5 km of the summit. One third of the interpolated volume is deposited beyond the 3 cm isopach line, highlighting the need to better constrain the thinning of the distal deposit. The volume of the cone formed with the N crater of OL (3–4 10⁶ m³; Kervyn et al. 2010) is not accounted for in this interpolation and should be added. This returns a total estimated volume of ~2. 10⁷ m³ for the 2007–08 eruption of OL, corresponding to an eruption magnitude of 3.5 (Pyle 2000).

Fig. 10

Interpolated ash layer thickness for the 2007–08 Oldoinyo Lengai eruption. Ash thickness is interpolated using an inverse distance weighting algorithm using thickness data collected in the field (see Fig. 9e), an isopach line of 3 cm around the ash-affected vegetation as obtained from the Principal Component Analysis (Fig. 9) and a constrain of extent of the 1 mm isopach line. Source of background: ArcGIS Map Services

Discussion

Potential and limitations of the mapping method

The method developed in this paper to map the impact of ash fallout on vegetation and vegetation recovery is based on methods to determine the impacts of forest fire (Diaz-Delgado et al. 1998; Salvador and Pons 1996; Riaño et al. 2002; Diaz-Delgado et al. 2003; Gouveia et al. 2010; Hope et al. 2012). Although the disruption process is different, the impact on vegetation is comparable yet more variable: according to the ash thickness plants might be totally buried, killed but with a rapid plant recovery due to seeds in the soil still being able to germinate (Anderson et al. 2012), or only covered temporarily by ash with a minor decay for a short time period (Ayiris and Delmelle 2012). In the latter case, the effect would be too limited to be detected by temporal change analysis of vegetation by remote sensing.

Although the zone affected by significant ash fallout is restricted to an ~100 km² area with low population density, the impact of the ash in terms of population relocation was more widespread. A survey of impacts on the population highlighted that thin ash fallout, which was insufficient to have a long term impact on the vegetation status, was enough to prevent livestock from grazing and forced temporary resettlements.

In order to map the effect of ash fallout on vegetation, it is required to isolate the effect of the ash relative to other environmental factors that might affect the state of the vegetation. In a region with very homogeneous vegetation, it is possible to compare the affected vegetation with undisturbed vegetation in the vicinity, in order to remove the effect of other causes of observed variation. In OL area, the topographic and vegetation gradients are substantial. It was therefore required to consider a long enough time series to understand the natural variation of the vegetation: as the variability of rainfall is the main factor controlling the state of vegetation (Reed et al. 2009), it is useful to remove this signal from the time series before mapping the vegetation change due to ash fallout. The limited knowledge of the rainfall pattern in the region led to the implementation of a simple rainfall correction method based on a single interpolated rainfall time series and specific vegetation response time for different classes of vegetation. Although this model is not ideal, it enables to partially remove the rainfall signal in the NDVI time series, and to better highlight the ash fallout pattern. Both the bi-temporal and PCA methods were able to map the extent of the ash affected zone. The bi-temporal comparison was most effective when applied at the peak of the growth season, when the contrast between affected and preserved vegetation was strongest. PCA analysis enabled to isolate the impact of ash on vegetation mainly in one component.

Recovery analysis was carried out on NDVI not corrected for rainfall variation, but these analyses were systematically implemented over time series of several years, in order to limit the influence of inter-annual seasonal variations. Both the trend analysis and the different recovery indices analysed had the ability to map the dynamic of vegetation change just after and within the five years following the end of the eruption. Results obtained at the end of the eruption highlight that rainfall-correction is not essential to identify ash-affected area, although the contrast with non-affected areas is enhanced thanks to this rainfall correction.

The difference between the different recovery methods lies in the contrast between the recorded evolution within the ash-affected zone and the variability caused by other factors in the rest of the study area. In that perspective the Change of NDVI (Tc) index proposed by Lu et al. (2012) performed best as it accounts for the normal seasonal variation in the time series.

The affected vegetation around OL is mostly characterized by a grassland and shrub savannah with scattered trees (Reed et al. 2009). Forests being only present on the flank of the Highlands were not significantly affected by the ash fallout. Such low vegetation is more likely to be affected by a limited ash thickness. Wilson et al. (2011) noted that in Argentina, 10 cm of ash was sufficient to totally bury winter pastures, but that vegetation could survive or quickly re-establish through an ash layer of <7 cm – or even <12 cm were there was no water shortage. Water and wind help vegetation recovery by stimulating erosion (Wilson et al. 2011; De Rose et al. 2011). Anderson et al. (2012) similarly highlighted the role of rainfall amount in controlling the germination of seeds after fires in the Serengeti region.

Our comparison with measured ash thickness suggests that 3 cm of ash was enough to disturb the state of the vegetation around OL sufficiently to be evidenced in our change analysis. However, zones affected by a thin ash layer were shown to recover over one year, whereas locations with more than 8–10 cm of ash needed several years to start recovering. It remains to be tested if a similar vegetation change analysis would perform well for forest vegetation. Similarly, the OL eruption took place in a dry environment with limited rainfall during the eruption: the role of rainfall in removing ash and therefore reducing impact on vegetation should also be considered when implementing this method in other climatic settings.

Vegetation recovery as a proxy for ash distribution?

Wilson et al. (2011) found a positive linear relationship between ash depth and length of land abandonment following an eruption. It was thus here postulated that the duration of the vegetation recovery or the recovery rate could be related to the ash thickness deposited on the vegetation. The elliptical shape and concentric pattern for the slope of the logarithmic trend of the recovery index over time and of the recovery duration indeed suggest a systematic variation of these parameters with ash thickness (Fig. 9). A power law statistical relationship is shown to approximate, although with low level of fitting, the relationship between measured ash thickness in the field and the slope value of the recovery trend. This suggests that the recovery rate is partially controlled by the ash thickness but that this relationship is non-linear and may be controlled by threshold values.

The limited spatial extent of the field measurements of ash thickness and the 250 m spatial resolution of the vegetation indices limits the accuracy of this statistical analysis and the identification of an ash thickness threshold sufficient to significantly affect the vegetation. Several centimeters of ash (3–5 cm) are probably enough to kill the surface vegetation, but seeds and roots are preserved in the soil and can lead to vegetation establishment quickly after the end of the eruption. On the other hand, for thicker ash coverage, the formation of a new soil layer is required to establish new plants, leading to a much longer recovery period.

As noted by previous authors (De Rose et al. 2011), vegetation recovery also depends on terrain morphology, geology and watershed morphology: as topography and water availability are characterized by large gradients in the study area, it can be expected that ash thickness will only be one factor controlling the recovery rate. Notwithstanding these limitations, the link between ash thickness and vegetation recovery enables to further constrain the variation in thickness of the ash layer around OL. Using conservative assumptions on the maximal extent of the ash layer, which cannot be inferred from vegetation change analysis, the volume derived in this study is two times larger than the 10⁷ m³ volume estimated by Kervyn et al. (2010) based on limited field reports. Although the distal layer is poorly constrained, the field observations constrain the ash layer on OL lower flank and base, and the vegetation analysis extends this constraint to the ~100 km² most affected by the ash fallout, including 64 % of the total estimated ash volume. The produced isopach map is an approximation of the ash deposited over the entire eruption period, and does not account for individual ash clouds which sometimes drifted in contrasted directions, depositing thin ash layers (i.e., <1 cm) to the north or east of the volcano.

Conclusion

We presented a temporal analysis of an 8-year long NDVI time series in the region around Oldoinyo Lengai volcano spanning across the 2007–08 eruption period. After correcting the time series for errors and rainfall-induced variation, a ~100 km² area where vegetation was severely impacted by ash was highlighted using bi-temporal comparison and principal component analysis. Linear temporal trend analysis and different vegetation recovery indices inspired from the literature on forest fire recovery were used to document the spatio-temporal pattern of vegetation recovery. The recovery rate and total duration was shown to be influenced by the spatial variation of ash thickness, with a thickness of 3 cm of ash being identified as the threshold for causing severe and systematic decay and subsequent recovery of the vegetation cover.

The proposed approach to estimate ash layer distribution and approximate thickness based on impact on the vegetation bears promises for characterizing ash layers in remote locations where field work is limited. Sufficient field observations are however required to calibrate the relationship. It is also expected that such an approach will work best for grassland and shrubland vegetation and in cases where rainfall did not remobilize the ash layer during or directly after the eruption.

The management of the OL eruption crisis and its impact by the Tanzanian authorities was limited to a restriction of access to the volcano during the eruption, without specific enforcement measures. No systematic effort was done in order to monitor the impact of the ash on the population and the environment, nor were recommendations or logistic support provided for the relocated population, to the exception of NGO’s. The results of vegetation change analysis presented here, if started directly after the eruption and continued over several months to years, would help to estimate the extent and duration of the impact and to direct long term recovery support to the communities located in the most affected zones as identified by the change detection analysis.