A Phenology-Based Method for Monitoring Woody and Herbaceous Vegetation in Mediterranean Forests from NDVI Time Series

Abstract

We present an efficient method for monitoring woody (i.e., evergreen) and herbaceous (i.e., ephemeral) vegetation in Mediterranean forests at a sub pixel scale from Normalized Difference Vegetation Index (NDVI) time series derived from the Moderate Resolution Imaging Spectroradiometer (MODIS). The method is based on the distinct development periods of those vegetation components. In the dry season, herbaceous vegetation is absent or completely dry in Mediterranean forests. Thus the mean NDVI in the dry season was attributed to the woody vegetation (NDVI_W). A constant NDVI value was assumed for soil background during this period. In the wet season, changes in NDVI were attributed to the development of ephemeral herbaceous vegetation in the forest floor and its maximum value to the peak green cover (NDVI_H). NDVI_W and NDVI_H agreed well with field estimates of leaf area index and fraction of vegetation cover in two differently structured Mediterranean forests. To further assess the method’s assumptions, understory NDVI was retrieved form MODIS Bidirectional Reflectance Distribution Function (BRDF) data and compared with NDVI_H. After calibration, leaf area index and woody and herbaceous vegetation covers were assessed for those forests. Applicability for pre- and post-fire monitoring is presented as a potential use of this method for forest management in Mediterranean-climate regions.

1. Introduction

Monitoring the dynamics of different vegetation components of the forest is essential to identify and understand trends in vegetation structure and ecosystem functioning. Specifically, changes in vegetation cover from woody to herbaceous species might indicate successional processes [1] or decline/recovery processes following natural and/or anthropogenic disturbances [2,3]. Woody and herbaceous vegetation constitute different layers in forests (e.g., under and over-story layers), which makes their monitoring important for use in land-atmosphere models [4]. Moreover, practical fields such as the management of wood [5] and range [6] production, fire risk [7], wildlife habitat [8] and biodiversity [9] depend upon estimating woody and herbaceous vegetation in woodlands. For example, woody and herbaceous vegetation regulate fire behavior and fire severity in different manners [10]. They also differ in their value for livestock grazing depending on animal type (e.g., cattle and sheep consume mainly herbaceous plants while goats feed mostly on woody species).

Conventional monitoring through ground inventory surveys is costly, time consuming and limited in temporal and spatial coverage [11]. In particular, reliable long-term vegetation monitoring in Mediterranean forests is challenging, as these forests are usually unable to return their costly management and monitoring operations. Satellite remote sensing offers an efficient and accessible tool for continuous vegetation monitoring at local to global scales complementing spatio-temporal limitations of field data.

Among satellite remote sensing tools, spectral vegetation indices are widely used for monitoring vegetation [12]. However, satellite-derived vegetation indices are mostly used to interpret changes such as greening or browning of the entire vegetation system without being able to distinguish between the woody and herbaceous contribution to those trends [13,14,15,16]. To date, a few studies have used vegetation indices to detect vegetation components in woodland systems at a sub-pixel level.

Canisius and Chen [17] proposed a physical-based approach to distinguish between understory and overstory reflectance in forests using satellite images. They applied multi-angle observations of canopy reflectance to monitor changes in understory reflectance. Yang et al. [18] proposed a more empirical-based method using the Moderate Resolution Imaging Spectroradiometer (MODIS) Bidirectional Reflectance Distribution Function (BRDF) product to retrieve ecosystem understory Normalized Difference Vegetation Index (NDVIu). Although both methods successfully retrieve understory reflectance in sparse forests [19] they are unable to separate the contribution of woody and herbaceous vegetation to the understory signal.

To distinguish between woody and herbaceous contributions, Roderick et al. [20] used time series analysis on NDVI. They aimed to map woody and herbaceous cover across Australia by applying a moving average method to separate between the long inter-annual trend and the seasonal (intra-annual) components of the time series. They assumed that woody vegetation in Australia is mainly evergreen (i.e., woody plants that keep their leaves and stay green throughout the entire year) and herbaceous vegetation is ephemeral (i.e., above ground tissues that wither during the dry season). Then, they attributed the baseline of the trend to the woody vegetation and the seasonal signal to herbaceous vegetation in the understory [20].

Lu et al. [21] further refined their method by applying a Seasonal-Trend decomposition (STL) [22] based on LOcally wEighted regreSsion Smooth (LOESS), adding a small seasonal component to the trend. By applying STL on 8-km NDVI time series derived from the Advanced Very High Resolution Radiometer (AVHRR), they produced ensemble-monthly averages of herbaceous vegetation cover for a 14-year period over Australia. Although promising, their proposed model is relatively complex and not suited for Mediterranean vegetation systems. In their model, the NDVI attributed to the woody vegetation increases simultaneously with the NDVI attributed to the herbaceous vegetation, while in actuality there is a delay of a few months between the development of ephemeral (herbaceous) and evergreen (woody) foliage in Mediterranean systems [23].

Here we present an efficient method to distinguish between NDVI from woody (i.e., evergreen) and herbaceous (i.e., ephemeral) vegetation in Mediterranean woodlands at a sub-pixel scale using phenology-based time series analysis. We first built 14-year MODIS NDVI time series at a temporal resolution of 16 days and spatial resolution of 250 m. Then, using the distinct phenology of the herbaceous and woody vegetation we decomposed NDVI time series into their woody (NDVI_W) and herbaceous (NDVI_H) signals. Results were compared with field data of overstory (i.e., woody) leaf area index and estimated woody and herbaceous vegetation covers in two Mediterranean woodlands. We also compared our results with NDVIu retrieved from MODIS BRDF data [18] at those two sites. Finally, we demonstrate the potential applications of this method for forest management and pre/post-fire monitoring.

2. Data and Methods

2.1. Satellite Data and Processing

We used the MODIS NDVI product (MOD13Q1) at a spatial resolution of 250 m. This product is a composite based on a single day value selected from periods of 16 days using a maximum value criterion [24]. Fourteen-year time series (from September 2000 to September 2014) of 16-day NDVI were built and woody and herbaceous signals were separated using the procedure described in Section 2.3 for each pixel at two woodland sites (see Section 2.2). NDVI, the most used vegetation index, is based on the solar reflectance in the red (Red, 0.6 μm) and near infrared bands (NIR, 0.8 μm). It is normalized to get values between −1 to 1, with fully vegetation cover approaching 1 [25]:

$$\mathrm{NDVI} = \frac{\mathrm{NIR}-\mathrm{Red}}{\mathrm{NIR} + \mathrm{Red}}$$

(1)

NDVI exhibits asymptotic (saturation) problems over highly dense vegetation cover and is also sensitive to canopy background variations with brighter canopies negatively biasing NDVI values over sparse vegetation [24]. However, NDVI reliably detected green biomass and vegetation dynamics in Mediterranean ecosystems that display relatively low to moderate values [26,27,28,29]. We preferred NDVI to the enhanced vegetation index (EVI), although EVI overcomes some of the above-mentioned drawbacks, because EVI tends to be much noisier over our study area with lower signal to noise ratio.

Noise and uncertainties in NDVI time series were reduced with LOESS [30], a technique for eliminating outliers from datasets with seasonal signal [21]. LOESS was applied with a window of n = 16 time steps (×16 days) to track the intra-annual change eliminating outliers that were greater than the standard deviation. Then, the residual was added to the smoothed time series to reproduce the original time series excluding only erroneous data. This procedure successfully eliminated cloud-contaminated and other erroneous data from the time series in all pixels within the study sites area.

To retrieve understory NDVI, we applied the method proposed by Yang et al. [18] on the MODIS BRDF/Albedo product (MCD43A1, version 5). The MCD43A1 is produced every eight days with 16 days of acquisition using data from both Terra and Aqua satellites [31]. It provides the weighting parameters associated with the RossThick–LiSparse BRDF model that describes the reflectance anisotropy at a 500m spatial resolution. The MCD43A1 product for tile h20v05 was downloaded with the associated data quality (MCD43A2) for the period of August 2005–August 2006.

The rationale behind Yang et al. [18] method is that pixels with smaller overstory LAI (LAIo) expose more understory cover making the understory contribution to the NDVI more significant. Thus NDVI in pixels with LAIo approaching 0 would be attributed mostly to the understory. Such pixels are difficult to find or even identify in forests. However, if LAIo varies between pixels, each pixel exposes different understory cover displaying distinct NDVI value. Then, the expected NDVI for LAIo = 0 (i.e., NDVIu) could be calculated by using regression analysis on those pixels and their reconstructed NDVI values for different angles (from BRDF data). A linear regression is performed after setting the nadir-view NDVI as the independent variable and the reconstructed NDVI values from the bidirectional angles as dependent variables. NDVIu is then calculated as the mean NDVI value in which regression lines (one regression line for each angle in the pixels of a defined window) have the smallest standard deviation. Correlations should be significant for pixels within a defined window with a similar canopy structure but variable LAIo. For the full description of this method and some illustrative examples the reader is referred to [18].

Here, we chose a window of 3 × 3 pixels following the criteria established by Yang et al. (see Section 5.2 in [18]). A 16-day NDVIu time series was produced for each one of the two case study sites. It is worth noting that NDVIu includes the combined signal from evergreen woody vegetation and ephemeral herbaceous plants in the understory. Therefore, NDVIu should be comparable to soil background NDVI (0.1–0.2) in Yatir during the dry season. This is because there is hardly any evergreen vegetation in the understory and the herbaceous vegetation is completely dry (see Section 2.2.1). Hence, comparable NDVIu and ecosystem NDVI values during the wet season in Yatir would imply that the herbaceous vegetation is the main contributor to the NDVI during this period assessing the method’s assumptions (see Section 2.3). At Mt. Carmel, NDVIu is expected to be higher than soil NDVI during the dry season due to evergreen woody species in the understory (see Section 2.2.2.).

2.2. The Case Study Sites

2.2.1. Yatir Pine Forest

The Yatir forest is a pine forest (Pinus halepensis) planted mostly from 1964 to 1969 in the semiarid region of Israel (31°20′49.2″N, 35°03′07.2″E, Figure 1a). It covers an area of ~2800 ha (Figure 1b). The average elevation is 650 m a.s.l. and mean annual precipitation is 285 mm∙yr⁻¹ (for the last 40 years). The mean annual temperature is 18.2 °C with 13 and 31 °C for mean winter (November–January) and summer (May–July) temperatures, respectively. Tree density is ~300 trees∙ha⁻¹ with a mean tree diameter of 17.5 cm and mean tree height of 11 m [32]. The canopy leaf area index (LAI) is typically ~1.4 ± 0.4 m²∙m⁻² with small fluctuations between winter and summer [33]. The understory is comprised of ephemeral herbaceous species (i.e., theropytes, geophytes and hemicryptophytes) growing during the wet season (September–April) and drying out in the beginning of the dry season (May–June). A relatively thin needle litter layer covers the forest floor during the needle senescence period (June–August) [23]. The forest understory is absent of evergreen vegetation (Figure 1c).

Figure 1. (a) Location of the two study sites (Yatir forest and Mt. Carmel woodlands). (b) Aerial (Google Earth ^®) and (c) field views of the planted pine forest of Yatir in the semiarid region of Israel. Photo: Eugene Ivanov.

Figure 1. (a) Location of the two study sites (Yatir forest and Mt. Carmel woodlands). (b) Aerial (Google Earth ^®) and (c) field views of the planted pine forest of Yatir in the semiarid region of Israel. Photo: Eugene Ivanov.

This forest site was selected because it has only two vertical layers (i.e., evergreen trees in the overstory and ephemeral herbaceous species in the understory) and a monospecific single-aged stand structure, which makes easier interpretation of the results. Moreover, this site was subject to extensive studies since the operation of its FLUXNET station in 2000 and is of special interest for its contribution, along other semi-arid ecosystems, to the climate system [34,35].

We used allometric mean annual LAIo estimates for the period of January 2000–June 2005 and monthly effective LAIo measured using a TRAC device (Tracing Radiation and Architecture of Canopies) during 2005 and 2006, both from Sprintsin et al. [33], to calibrate and compare with NDVI_W (i.e., overstory trees in the case of Yatir).

2.2.2. Mt. Carmel Mixed Pine-Oak Woodlands

Mt. Carmel region (35°E, 32°N; 10–520 m, Figure 1a) has a typical East-Mediterranean climate with dry hot summers (June–August) and cool rainy autumns and winters (October–May). The mixed pine-oak woodlands in Mt. Carmel cover an area of ~21,500 ha (Figure 2a). Mean annual precipitation ranges from 550 mm∙yr⁻¹ near the coastal plane to 750 mm∙yr⁻¹ at the highest elevations (550 m a.s.l.). Vegetation in Mt. Carmel varies between dense multi-aged mixed pine-oak woodland and a more open and patchy shrublands. It consists of a variety of woody shrub and tree species in which the most dominant are Pinus halepensis, Quercus calliprinos, Pistacia lentiscus L. and Cistus salviifolius L., along with a large variety of herbaceous species (i.e., therophytes, geophytes and hemicryptophytes). This results in a complex vertical and horizontal multi-strata canopy structure [36] (Figure 2b).

Figure 2. (a) Aerial (Google Earth ^®) and (b) field views of the mixed pine-oak evergreen woodlands at Mt. Carmel. The burnt area from the wildfire of 2010 (red line) and the location of the 22 survey plots (dots) are indicated in (a). Photo: Naama Tessler.

Figure 2. (a) Aerial (Google Earth ^®) and (b) field views of the mixed pine-oak evergreen woodlands at Mt. Carmel. The burnt area from the wildfire of 2010 (red line) and the location of the 22 survey plots (dots) are indicated in (a). Photo: Naama Tessler.

The most dominant tree species are the Aleppo Pine (P. halepensis) and the common oak (Quercus calliprinos). These two evergreen tree species vary strongly in their total cover and shift in dominance depending on ecological site conditions (mainly bedrock type), land use history (grazing, cutting), and disturbances (fire). Herbaceous plants in the understory are all ephemerals with a life cycle during the wet season. The area went through repeated fires during the last 20 years affecting the regeneration of P. halepensis, which resulted in multi-aged pine distribution [37]. Wildfires in this site changed the forest structure and ratio of woody to herbaceous vegetation cover affecting its hydrological balance and ecosystem functioning [38].

We selected this site to test our method in a more complex woodland system comprising evergreen woody vegetation in the understory. In addition, this area went through a large wildfire (~2500 ha) in December 2010 [39] and was monitored since then providing in situ data of vegetation cover that can be used to compare with our results. This site represents a typical complex-structured Mediterranean woodland system with relatively high productivity (unpublished data).

Four field surveys were conducted since the wildfire of 2010 during the spring seasons (March–May, i.e., the time of peak herbaceous cover) of 2011–2014 in 22 plots of 10 × 10 m² each. The plots were selected based on Mt. Carmel geology map on southern and northern aspects within the burnt area (Figure 2a). In each plot, total vegetation cover and specific species cover were determined following Stohlgren [40]. Species cover were summed and grouped into three groups: evergreen trees in the overstory (Tr-Ev), evergreen shrubs in the understory (Sh-Ev) and ephemeral herbaceous plants in the understory (He-Ep) (Table 1). In addition to vegetation cover assessment, fire severity was determined in each plot and classified as: low, medium and high according to the classification proposed by Neary et al. [41]. Classification was extended to the total burnt area by using high spatial resolution aerial photos and GIS tools.

Table 1. The mean vegetation cover (%) of evergreen trees (Tr-Ev) and shrubs (Sh-Ev) and ephemeral herbaceous species (He-Ep) estimated in 22 plots at Mt. Carmel (Figure 2a) during 2011–2014. The total vegetation cover in each plot, and the mean and standard errors (SE) of all plots are also indicated.

Plot	Tr-Ev (over)	Sh-Ev (under)	He-Ep (under)	Total (over + under)
1	14	37	17	68
2	9	38	33	80
3	2	39	44	85
4	3	25	59	87
5	12	41	30	83
6	28	41	16	85
7	15	30	30	77
8	17	46	18	82
9	16	29	23	68
10	16	34	20	70
11	9	53	26	88
12	19	28	30	77
13	20	31	23	73
14	9	17	35	60
15	14	29	26	68
16	25	38	20	83
17	31	26	36	93
18	13	24	29	67
19	28	25	15	68
20	14	49	18	80
21	29	18	31	78
22	28	46	11	85
AVG	17	34	27	78
SE	2	2	2	2

Table 1. The mean vegetation cover (%) of evergreen trees (Tr-Ev) and shrubs (Sh-Ev) and ephemeral herbaceous species (He-Ep) estimated in 22 plots at Mt. Carmel (Figure 2a) during 2011–2014. The total vegetation cover in each plot, and the mean and standard errors (SE) of all plots are also indicated.

Plot	Tr-Ev (over)	Sh-Ev (under)	He-Ep (under)	Total (over + under)
1	14	37	17	68
2	9	38	33	80
3	2	39	44	85
4	3	25	59	87
5	12	41	30	83
6	28	41	16	85
7	15	30	30	77
8	17	46	18	82
9	16	29	23	68
10	16	34	20	70
11	9	53	26	88
12	19	28	30	77
13	20	31	23	73
14	9	17	35	60
15	14	29	26	68
16	25	38	20	83
17	31	26	36	93
18	13	24	29	67
19	28	25	15	68
20	14	49	18	80
21	29	18	31	78
22	28	46	11	85
AVG	17	34	27	78
SE	2	2	2	2

Average in situ vegetation cover for each plot were used to calibrate NDVI_W and NDVI_H, and to produce woody and herbaceous cover maps for 2000–2014 at Mt. Carmel. NDVI_H was regressed against He-Ep, while NDVI_W was regressed against the sum of Tr-Ev and Sh-Ev covers. Calibrations included only 14 MODIS pixels (250 m) because some of the 22 plots were within the area of a single pixel. In such case, species covers were averaged to get the mean cover within the pixel area. Plots that fall between pixels were discarded.

Although there are some uncertainties when comparing field and satellite data due to scaling and/or sensor view issues [42], this method has been extensively used to calibrate/validate satellite vegetation indices with biophysical parameters [21,26,27,43]. To test the robustness of this calibration we calculated vegetation covers using the classical two-end members fractional vegetation cover equation [44]:

$$\mathrm{FVC} = \frac{{\mathrm{NDVI}}_X-{\mathrm{NDVI}}_0}{{\mathrm{NDVI}}_{\mathrm{FULL}}-{\mathrm{NDVI}}_0}$$

(2)

where NDVI₀ in Equation (2) is the NDVI from the soil background (0.1), NDVI_Full is the NDVI from a fully vegetated area (0.7 and 0.9 for woody and herbaceous, respectively) and NDVI_X is the NDVI_W or the NDVI_H value. Finally, we compared woody and herbaceous cover maps retrieved from the two models.

2.3. Description of the Method

The well-defined wet and dry seasons characterizing Mediterranean climate regulate the photosynthetic activity and growth of the vegetation in Mediterranean ecosystems [45]. The seasonal year in Mediterranean forests is defined as the 12-month period since the beginning of the rainfall in October [32]. Ephemeral herbaceous species appear as an understory layer in forests shortly after the beginning of the rainy season (October–November) reaching a peak biomass at the end of winter (February), drying out in spring (April) [46]. In contrast, the evergreen woody vegetation becomes most active from early spring (March) developing new leaves towards the dry season (June–August) [23] (Figure 3).

Figure 3. Schematic representation of the distinct growth and senescence periods of evergreen woody vegetation (dashed blue line) and ephemeral herbaceous plants (red line) in Mediterranean forests. Phenological stages are shown as the relative Normalized Difference Vegetation Index (i.e., NDVI/NDVI_max) of each of those two vegetation components. The wet and dry periods are also indicated.

Figure 3. Schematic representation of the distinct growth and senescence periods of evergreen woody vegetation (dashed blue line) and ephemeral herbaceous plants (red line) in Mediterranean forests. Phenological stages are shown as the relative Normalized Difference Vegetation Index (i.e., NDVI/NDVI_max) of each of those two vegetation components. The wet and dry periods are also indicated.

These two components, i.e., woody and herbaceous species, contribute differently to the NDVI (Figure 3) resulting in a mixed signal at the sub pixel scale [19,20]. Hence, only by using time series analysis the contribution of these components to the total ecosystem NDVI (NDVI_Ecos) can be estimated. Because in Mediterranean woodlands the herbaceous vegetation is mostly ephemeral, NDVI_Ecos during the dry season (June–August) could be attributed solely to evergreen woody species and soil background (hereafter NDVI_W). Then, the seasonal component (i.e., the NDVI_Ecos in wet season—October–April) would be attributed mostly to the herbaceous vegetation (hereafter NDVI_Seas), which is the most dynamic component (in terms of NDVI) in those systems [21]. Accordingly, we separated the woody and herbaceous NDVI from the NDVI_Ecos in three steps:

(1)

The average NDVI_Ecos over the dry period (June–August) was calculated and taken as NDVI_W for each seasonal year (i.e., from September to September). If NDVI_Ecos in the wet season was lower than the calculated NDVI_W, the minimum NDVI_Ecos value was taken as the NDVI_W instead. This ensures that abrupt intra-annual changes following disturbances (e.g., fires or clearing) are detected. For example, if a fire event occurs in December (i.e., during the wet season), taking the average NDVI_Ecos over the dry season (June–August) would overestimate woody cover in that specific year. In contrast, taking the minimum NDVI_Ecos, which is the NDVI_Ecos value following the fire, would be more representative because it includes the change due to the fire.

(2)

The NDVI_W is then subtracted from NDVI_Ecos to compute the seasonal component of the time series (NDVI_Seas).

(3)

The maximum NDVI_Seas value in each seasonal year is taken as NDVI_H, which represents the peak biomass/green-cover of the herbaceous vegetation [13].

By assuming a constant NDVI for soil background during the dry season [27], NDVI_W could be used to track inter-annual changes in forest’s woody cover or leaf area index (LAI) after calibrating against ground truth data. Similarly, NDVI_H could be used to monitor changes in green biomass [26,27] or herbaceous vegetation fraction cover [20,21].

Although NDVI_H does not account for herbaceous vegetation fully covered by the canopy, it is assumed that in closed forests it contributes a negligible fraction due to the limited light resource [21]. In open forests, though, because sensor view angle changes throughout the year, NDVI_H is assumed to represent almost all of the herbaceous vegetation present in the forest floor. The validity of these assumptions were previously demonstrated in several studies [21,26,27] and will be further examined here by comparing with field data and the calculated understory NDVI (NDVIu) from the MODIS BRDF data [18].

2.4. Applications of NDVI_W and NDVI_H for Pre- and Post-Fire Monitoring in Mt. Carmel

We demonstrate three pre/post fire-monitoring applications using NDVI_W and NDVI_H in Mt. Carmel woodlands selecting the wildfire of 2010 as a case study. Field data collected at the burnt area (described in Section 2.2.2) were compared with three applications of NDVI_W and NDVI_H:

(i)

Post-fire monitoring of woody and herbaceous recovery (i.e., changes in woody and herbaceous cover) in the burnt area of the wildfire of 2010 (four post-fire years) assessed from NDVI_W and NDVI_H.

(ii)

Fire severity assessment in the burnt area from the difference in NDVI_W between pre and post fire years.

(iii)

Production of a fuel-based fire risk map from NDVI_W for the year prior to the wildfire and comparison with fire-spread behavior (i.e., the burnt area). Fire risk map was produced by assigning a relative score from 1 to 10 (i.e., from the lowest to the highest risk level) to each pixel in the Mt. Carmel area according to its woody cover (i.e., 1 for minimum and 10 for maximum cover) and dryness status (i.e., 1 for the most positive, or no trend, and 10 for the most negative trend).

Although forest fires are driven also by weather conditions, topography and many other factors, here we examined the importance of fuel density and status on fire spread behavior. This is important because fire hazard models mostly use static fuel-based maps that does not account for long-term drought effects on the woody vegetation status. We did not use NDVI_H in this map because fires in Mt. Carmel were mostly driven by woody dry matter and less by ephemeral herbaceous vegetation [47]. Moreover, the fire risk map produced when using also NDVI_H was mostly similar to that produced from NDVI_W alone.

3. Results and Discussion

3.1. Mapping Leaf Area Index from NDVI_W in Yatir

Figure 4a shows the time series of NDVI of the entire ecosystem (NDVI_Ecos) from 2000 to 2014 (from September to September) in one representative pixel at the planted pine-forest of Yatir. The NDVI_Ecos during the summer was relatively low (~0.3) increasing in winter, reaching a maximum value of 0.45–0.6 and decreasing afterward towards the next summer (Figure 4a). NDVI_W generally increased from 2001 (i.e., the seasonal year of 2000/1) to 2008 decreasing afterward till the end of the time series (Figure 4c).

NDVI_Seas was very dynamic throughout the years in Yatir (Figure 4e) probably due to variations in annual rainfall amount [26]. The NDVIu derived from BRDF data was almost the same as the NDVI_Ecos during the main wet season (November–April, Figure 5a). Because understory vegetation in Yatir is comprised mostly of ephemeral herbaceous species, this supports the assumption that herbaceous vegetation in the understory is the main contributor to the seasonal change in NDVI_Ecos. However, NDVI_Ecos increased in October before NDVIu started to increase, probably due to tree leaf development as observed from the increase in LAIo (Figure 5a). Similarly, NDVI_Ecos was higher than NDVIu during the senescence period (March–April) coincident with an increase in LAIo (Figure 5a).

The relatively high NDVIu during the summer (~0.23) compared to the soil background NDVI for this region (0.1–0.2) was probably due to accumulated greenish to yellowish needle litter in the understory, which is typically accumulated during this period in the forest floor [23]. NDVI_H showed almost an opposite trend to NDVI_W mostly decreasing till 2009, increasing in the following years (Figure 4g). The relationship between LAIo from the TRAC measurements and June to August NDVIE_cos was significant (R² = 0.72, p < 0.05, Figure 5b). NDVI_W also agreed well with in situ mean annual LAIo during January 2000–June 2005 (R² = 0.71, p < 0.05), exhibiting a typical exponential relationship (LAI = 0.81⋅e^7⋅NDVI_P, Figure 5c). Indeed, the significant relationships between NDVI and in situ LAIo imply that the mean over the dry season NDVI_Ecos (i.e., NDVI_W) is a good estimator of the mean annual LAIo in this forest site.

Figure 4. Examples of the original and decomposed NDVI time series in one representative pixel (250 m) at the evergreen pine forest of Yatir (left column) and the pine-oak woodlands of Mt. Carmel (right column). Original and smoothed time series of NDVI_Ecos are shown in (a,b), NDVI_W in (c,d), NDVI_Seas in (e,f) and NDVI_H in (g,h) for Yatir and Mt. Carmel, respectively.

Figure 4. Examples of the original and decomposed NDVI time series in one representative pixel (250 m) at the evergreen pine forest of Yatir (left column) and the pine-oak woodlands of Mt. Carmel (right column). Original and smoothed time series of NDVI_Ecos are shown in (a,b), NDVI_W in (c,d), NDVI_Seas in (e,f) and NDVI_H in (g,h) for Yatir and Mt. Carmel, respectively.

Figure 5. (a) NDVI_Ecos, NDVIu (retrieved from Bidirectional Reflectance Distribution Function product), and monthly in situ overstory LAI (LAIo) in a 4 km² area at the Yatir pine forest. Scatterplots showing the relationships between (b) monthly LAIo vs. NDVI_Ecos and (c) mean annual LAIo vs. NDVI_W (2000–2006) for the same area.

Figure 5. (a) NDVI_Ecos, NDVIu (retrieved from Bidirectional Reflectance Distribution Function product), and monthly in situ overstory LAI (LAIo) in a 4 km² area at the Yatir pine forest. Scatterplots showing the relationships between (b) monthly LAIo vs. NDVI_Ecos and (c) mean annual LAIo vs. NDVI_W (2000–2006) for the same area.

Figure 6 shows the spatial distribution of the mean annual LAIo calculated from NDVI_W in Yatir at a 250 m spatial resolution with mostly negative trends for the period of 2000–2014.

Figure 6. Spatial distribution of (a) the mean annual and (b) trends of overstory LAI retrieved from NDVI_W at Yatir for 2000–2014. Significant trends are indicated in (a) as + (positive) and • (negative), while in (b) all significant trends are indicated as +.

Figure 6. Spatial distribution of (a) the mean annual and (b) trends of overstory LAI retrieved from NDVI_W at Yatir for 2000–2014. Significant trends are indicated in (a) as + (positive) and • (negative), while in (b) all significant trends are indicated as +.

3.2. Assessing Woody and Herbaceous Cover from NDVI_W and NDVI_H in Mt. Carmel

Figure 4b,d,f,h show the NDVI_Ecos and its decomposition into NDVI_W, NDVI_H and NDVI_Seas in one representative pixel at the mixed pine-oak woodlands of Mt. Carmel. The large wildfire of December 2010 (2500 ha) is noted by an abrupt decline in NDVI_W from 2010 to 2011 (Figure 4d). Another, much smaller wildfire event in 2005 (150 ha, [39]) is also noted by a decrease in NDVI_W from 2004 to 2005. As in Yatir, the NDVIu retrieved from BRDF data suggests that understory vegetation is the primary contributor to NDVI_Ecos during the wet season (November–April, Figure 7). However, unlike in Yatir, the NDVIu was high in Mt. Carmel woodlands during the dry season (~0.4) due to the relatively dense evergreen vegetation cover in the understory.

Figure 7. NDVI_Ecos and NDVIu (from BRDF) in a 4-km² area of the pine-oak woodlands at Mt. Carmel.

Figure 7. NDVI_Ecos and NDVIu (from BRDF) in a 4-km² area of the pine-oak woodlands at Mt. Carmel.

Woody and herbaceous vegetation covers estimated in field were significantly correlated with NDVI_W and NDVI_H, respectively (p < 0.01 for both, Figure 8). The coefficient of determination was moderate with R² of 0.58 for NDVI_W vs. woody cover and 0.69 for NDVI_H vs. herbaceous vegetation cover. Correlations were only moderate probably due to the scale differences between the field plots (10 × 10 m²) and the MODIS pixels (250 × 250 m²).

Figure 8. Scatterplots of woody (left) and herbaceous (right) vegetation covers (%) assessed in field vs. NDVI_W and NDVI_H in 14-MODIS pixels. Each dot in the graph is the four-year averaged vegetation cover within one MODIS pixel (see Section 2.2.2). For the herbaceous vegetation cover, an exponential function was fitted following the assumption that NDVI_H equals 0 in the dry season when herbaceous vegetation is absent (0% cover).

Figure 8. Scatterplots of woody (left) and herbaceous (right) vegetation covers (%) assessed in field vs. NDVI_W and NDVI_H in 14-MODIS pixels. Each dot in the graph is the four-year averaged vegetation cover within one MODIS pixel (see Section 2.2.2). For the herbaceous vegetation cover, an exponential function was fitted following the assumption that NDVI_H equals 0 in the dry season when herbaceous vegetation is absent (0% cover).

We used the regression functions from Figure 8 to derive maps of the mean annual woody and herbaceous vegetation covers at Mt. Carmel for the period of 2000–2014 (Figure 9a,d). Similar spatial patterns with slightly different magnitudes were obtained from NDVI_W and NDVI_H when using the two-end members FVC equation (Equation (2), Figure 9b,e).

Figure 9. The 14-year mean woody (a–c) and herbaceous (d–f) vegetation cover (%) and NDVI trends in Mt. Carmel. Vegetation cover was estimated from NDVI_W and NDVI_H using (a,d) NDVI-field regression functions (see Figure 8) and (b,e) the two-end members FVC equation (Equation (2)). The strong NDVI_W decline in (c) and contrast NDVI_H increase in (f) are the result of the 2010 wildfire (compare with wildfire area in Figure 2a).

Figure 9. The 14-year mean woody (a–c) and herbaceous (d–f) vegetation cover (%) and NDVI trends in Mt. Carmel. Vegetation cover was estimated from NDVI_W and NDVI_H using (a,d) NDVI-field regression functions (see Figure 8) and (b,e) the two-end members FVC equation (Equation (2)). The strong NDVI_W decline in (c) and contrast NDVI_H increase in (f) are the result of the 2010 wildfire (compare with wildfire area in Figure 2a).

However, the range in woody cover at Mt. Carmel was narrower when using the calibration function (30%–70%, Figure 9a) compared to the obtained from the two-end members function (15%–80%, Figure 9b). This was probably because of the relatively narrow range of woody covers estimated at field and used for the calibration (40%–60%, Figure 8). Omitting low and high covers from the calibration resulted in overestimation and underestimation of the woody cover at sparse/dense woody cover, respectively.

Oppositely, the range of herbaceous cover was wider when using the calibration function (20%–90%, Figure 9d) compared to the retrieved from the two-end members function (30%–65%, Figure 9e). This is probably because of the exponential function used for calibration, whereas the function used in two-end members model is linear. Therefore, the two-end members model underestimated the herbaceous vegetation cover at low and high values of vegetation covers.

NDVI_W and NDVI_H were useful in assessing woody and herbaceous vegetation dynamics relatively, although the described discrepancies. A good example is the effect of the 2010 wildfire on the woody and herbaceous trends observed in Figure 9c,f. A large decline in NDVI_W is clearly noted in the 14-year NDVI_W trends (~0.2 per decade), while NDVI_H slightly increased (<0.1 per decade) during the same period.

3.3. Pre and Post Fire Assessment Using NDVI_W and NDVI_H in Mt. Carmel

NDVI is an indicator of the vegetation dryness [48] and could be used to assess dry woody vegetation matter, which is the main fuel source in Mediterranean forest wildfires. A decline in NDVI_W would imply drying trends in woody vegetation and increasing dry matter for potential fires [47]. On the contrary, low woody cover would reduce fire vulnerability in woodlands because there is less fuel for fire. Ephemeral herbaceous vegetation plays only a minor role in Mt. Carmel wildfires because of its relatively low biomass.

Figure 10a shows that the fuel-based fire risk map based only on the mean annual and trends in NDVI_W reproduced the fire-spread behavior of the 2010 wildfire with great accuracy.

Figure 10. (a) A fuel-based fire risk map produced for the year 2009 (prior to the 2010 wildfire) from the woody vegetation cover (mean NDVI_W) and its dryness status (NDVI_W trends). Superimposed is the area of the wildfire; histograms of (b) the total number of pixels in Mt. Carmel with their respective risk levels and (c) the ratio between the number of pixels with a specific risk level in the burnt zone to that in the entire Mt. Carmel area (in %). The dashed line in (c) indicates the percent area of burnt zone from the entire Mt. Carmel area (i.e., 8.6%).

Figure 10. (a) A fuel-based fire risk map produced for the year 2009 (prior to the 2010 wildfire) from the woody vegetation cover (mean NDVI_W) and its dryness status (NDVI_W trends). Superimposed is the area of the wildfire; histograms of (b) the total number of pixels in Mt. Carmel with their respective risk levels and (c) the ratio between the number of pixels with a specific risk level in the burnt zone to that in the entire Mt. Carmel area (in %). The dashed line in (c) indicates the percent area of burnt zone from the entire Mt. Carmel area (i.e., 8.6%).

Most of the pixels within the wildfire area were ranked with the highest risk level, being 18% and 21% (for levels 9 and 10, respectively) from the total cells ranked with the same levels in Mt. Carmel (Figure 10c). This is more than the expected in the burnt zone, which comprises only ~9% of the entire Mt. Carmel area. These results, based only on fuel density and status (both from NDVI_W), do not consider topographic or weather conditions, which are important driving forces of fire spread behavior. However, such a fuel-based map might improve fire hazard models that mostly rely upon static fuel maps [47]. Those do not consider the effects of drought years on the woody (and/or herbaceous in other cases) vegetation status, which is an important factor determining fire spread behavior as shown here.

Figure 11. Maps of (a) low, medium and high severity burnt areas classified in field and extended with high-resolution aerial photograph; and (b) the difference between post- and pre-fire NDVI_W (ΔNDVI_W). (c) Box plot of mean, first and third quartiles (with respective standard deviations) of ΔNDVI_W in the low, medium and high severity areas mapped at field (shown in a). Different letters indicate statistically significant differences at p < 0.001 using a two-tailed Student’s t-test, after a Bonferroni correction.

Figure 11. Maps of (a) low, medium and high severity burnt areas classified in field and extended with high-resolution aerial photograph; and (b) the difference between post- and pre-fire NDVI_W (ΔNDVI_W). (c) Box plot of mean, first and third quartiles (with respective standard deviations) of ΔNDVI_W in the low, medium and high severity areas mapped at field (shown in a). Different letters indicate statistically significant differences at p < 0.001 using a two-tailed Student’s t-test, after a Bonferroni correction.

NDVI_W was also used to assess fire severity in the burnt area of Mt. Carmel. Figure 11 shows the difference between NDVI_W from the year prior (2010) and posterior (2011) to the wildfire (ΔNDVI_W, Figure 11b). The mean ΔNDVI_W was significantly different (p < 0.001) between areas classified in field as low, medium and high severity with mean ΔNDVI_W of 0.1, 0.13 and 0.16, respectively (Figure 11c).

Finally, we used NDVI_W and NDVI_H to assess woody and herbaceous vegetation recovery (i.e., changes in cover) in the burnt area (Figure 12). Results from NDVI_W and NDVI_H agreed well with the field estimates when using both, the empirical relationships (from Figure 8) and the classical two-end members equation (Equation (2)). All estimates showed that after a woody cover reduction due to the wildfire (Figure 9c and Figure 11b), woody vegetation cover decreased from 2011 to 2012 (from the year after the wildfire to the next year) but increased in the subsequent two years (Figure 12).

Figure 12. Changes in woody and herbaceous vegetation cover following the wildfire of 2010 as estimated from field (open bars) and from NDVI (solid bars) using a field-based calibration function (NDVI-field) and two-end members fraction of vegetation cover (FVC) equation. Asterisks indicate statistically significant differences between the NDVI’s components and the field estimates for each specific year at p < 0.05. Different letters denote statistically significant differences in vegetation covers between years at p < 0.05.

Figure 12. Changes in woody and herbaceous vegetation cover following the wildfire of 2010 as estimated from field (open bars) and from NDVI (solid bars) using a field-based calibration function (NDVI-field) and two-end members fraction of vegetation cover (FVC) equation. Asterisks indicate statistically significant differences between the NDVI’s components and the field estimates for each specific year at p < 0.05. Different letters denote statistically significant differences in vegetation covers between years at p < 0.05.

The herbaceous vegetation cover that generally increased in the burnt zone (Figure 9f), showed an opposite change to the woody cover. It first increased from 2011 to 2012 (first to second year after the wildfire) but decreased afterward during the next two years (in 2013 and 2014).

The agreement between field estimates and calculated vegetation covers from NDVI_W and NDVI_H demonstrate the potential use of this method to monitor succession processes in forests following disturbances.

4. Conclusions

In this study, we proposed a phenology-based method that uses the distinct timing of ephemeral herbaceous (i.e, above ground tissues that wither during the dry season) and woody (i.e., evergreen) foliage development to separate time series of the Normalized Difference Vegetation Index (NDVI) into their respective contributions. The average NDVI during the dry season (i.e., when herbaceous species are dry or absent) was attributed to the woody vegetation (NDVI_W), while the maximum value during the wet season (NDVI_H) was attributed to the herbaceous vegetation after subtracting NDVI_W from the time series. The assumption that the herbaceous vegetation contributes most of the NDVI signal during the wet season was validated with the understory NDVI (NDVIu) retrieved from an independent semi-empirical method based on the small bidirectional variation in NDVIu compared to the canopy-level variation [18]. Both methods showed that change in NDVI during the wet season were caused by the development of herbaceous vegetation at the understory.

NDVI_W agreed well with the in situ monthly and mean annual overstory LAI (LAIo) measured at the planted pine forest of Yatir (southern Israel). The coefficients of determination (R²) were 0.72 and 0.71 for monthly and annual LAIo, respectively. NDVI_W and NDVI_H were also moderately correlated to woody and herbaceous vegetation covers estimated at field in Mt. Carmel pine-oak woodlands (northern Israel) with R² of 0.58 and 0.69, respectively (p < 0.1 for both correlations). The relationships between NDVI_W and LAIo were used to assess spatial variations in mean annual and trends in LAIo at Yatir for 2000–2014. The relationships between NDVI_W/NDVI_H and in situ woody/herbaceous covers were used to produce maps of mean annual and trends in woody and herbaceous vegetation covers at Mt. Carmel for 2000–2014 (using the calibration functions). The spatial patterns of those maps were comparable to those retrieved using a simple two-end members equation.

To show the potential use of this method, we presented three pre/post-fire applications using the Mt. Carmel 2010 wildfire as a case study: (1) We produced a fuel-based fire risk map with mean NDVI_W representing woody cover and its trend woody dryness condition (i.e., negative trends meaning drying trends). Pixels with high NDVI_W and negative NDVI_W trends were ranked as high fire-risk pixels and vice-versa. The fire risk map reproduced fire-spread behavior of the wildfire of 2010 with great accuracy. About 39% of the pixels in Mt. Carmel with the highest risk levels were found within the area of the wildfire, which constitute only 9% from the entire Mt. Carmel area; (2) We used the difference in NDVI_W between the year after and before the wildfire (ΔNDVI_W) to map fire severity. The mean ΔNDVI_W was significantly different (p < 0.001) between areas classified at field as low, medium and high severity with, respectively, decreasing ΔNDVI_W values; (3) We used NDVI_W and NDVI_H to monitor the recovery of the woody and herbaceous vegetation in the burnt area during post-fire years (2011–2014).

The proposed method has the potentials to be used in Mediterranean systems for several forest management, fire risk assessment and post fire monitoring purposes. Moreover, it could be used to assess overstory and understory LAI and fraction of absorbed photosynthetic active radiation (fPAR) [49], which are essential parameters in land surface modeling. Yet, it should be noted that limitations of satellite imaging in detecting understory vegetation underneath a denser overstory cover (>80%) must be further determined. This could be done, for example, by placing ground spectral sensors at several heights in the forest to determine over- and under-story NDVI [50] and compare with NDVI_W and NDVI_H.

Once established, this method holds a great potential for monitoring, understanding and managing complex vegetation systems comprising various combinations of woody and herbaceous vegetation in Mediterranean ecosystems.