**2. Methodology**

356 International Perspectives on Global Environmental Change

Remote sensing tools can be used to detect both evident functional changes produced by land-use transformations (Volante et al., In press), and other subtle and less noticeable changes including insect outbreaks (Kharuk et al., 2003), wind (Yuan et al., 2002), droughts (Tucker & Choudhury, 1987) or floods (Sanyal & Lu, 2004), fires (Riano et al., 2002), pollution (Chu et al., 2003), etc. These impacts may derive in significant changes in key ecological processes, for instance, carbon balance, microclimate, and biodiversity patterns (Turner, 2005; Lovett et al., 2006; Perry & Millington, 2008). Remote sensing has been proved to be useful for monitoring this kind of "within-state" changes (Vogelmann et al., 2009). In particular, satellite-derived spectral vegetation indices, such as the Enhanced Vegetation Index (EVI) and the Normalized Difference Vegetation Index (NDVI), are considered the most useful approach to monitor ecosystem responses to environmental changes (Pettorelli et al., 2005). Vegetation indices constitute the most feasible approach to estimate primary production at the regional scale (Paruelo et al., 1997) since they show a linear response to the intercepted fraction of photosynthetically active radiation (FPAR) (Hanan et al., 1995), which represents the conceptual basis to relate vegetation indices with net primary

Where NPP is the Net Primary Production, PAR is the amount of incident Photosynthetically Active Radiation, FPAR is the fraction of that PAR that is intercepted by vegetation green tissues, and RUE is the Radiation Use Efficiency that plants have to transform that radiation into organic carbon compounds. Given this direct relationship with NPP, the most integrative descriptor of ecosystem functioning (McNaughton et al., 1989; Virginia & Wall, 2001), vegetation indices are frequently used to derive indicators of ecosystem functioning such as the annual amount of carbon absorbed by vegetation, or the seasonality and phenology of the carbon gain dynamics (Pettorelli et al., 2005; Alcaraz-

To evaluate the usefulness of satellite-derived vegetation indices for monitoring functional changes within protected areas, we focused on the Sub-Mediterranean Pyrenean oak forests (*Quercus pyrenaica* Willd.) of the Sierra Nevada National Park (Spain). These forests are considered as a Natural Habitat of Community Interest (*Quercus pyrenaica* oak woods and *Quercus robur* and *Quercus pyrenaica* oak woods from Iberian northwestern, Directive 92/43/CEE) (García & Mejías, 2009). The Pyrenean oak forests are a *quasi*-endemic habitat of the Iberian Peninsula. The only non-Iberian representations are in the Central West of France and in the Rif Mountains of northern Morocco. In the South of Spain, the Pyrenean oak is considered as a vulnerable species (Blanca & Mendoza, 2000). Sierra Nevada oak populations are considered of great biogeographical importance since they constitute the southernmost Iberian representation of these forests (Molero et al., 1992) and they are considered relict deciduous forests in the Southern Mediterranean region (Blanca & Mendoza, 2000; Blanca, 2001). Several stands of these forests in the Sierra Nevada National Park have an unfavorable conservation status (Molero et al., 1992; Bonet et al., 2010). Multiple global change drivers have an impact on these southernmost woodlands of *Quercus pyrenaica* in the Iberian Peninsula. Historically, these populations have been subjected to intense human disturbances (logging, fires, grazing, agriculture, etc). As a result, these forests are highly fragmented and display low ecological maturity (García & Mejías, 2009) that threatens their long-term conservation. Currently, trends towards temperature rises and precipitation decreases have been hypothesized as the main constraining factor reducing

NPP = PAR \* FPAR \* RUE (1)

production (NPP) through Monteith's model (Monteith, 1972) (equation 1).

Segura et al., 2006).

#### **2.1 The Pyrenean oak forests of Sierra Nevada National Park**

Sierra Nevada National Park is located in the southeast of the Iberian Peninsula (Figure 1). This National Park protects the best samples of high and medium Mediterranean mountainous ecosystems (MMARM, 2004). This park is a hot spot for plant species richness (Blanca et al., 1998; Blanca, 2001) and invertebrate biodiversity. Its altitude (several summits over 3000 m.a.s.l.), its proximity to Africa, and steep altitudinal gradient constitute the main ecological and evolutionary factors determining its high biodiversity.

The Pyrenean oak forests (Figure 1) of Sierra Nevada represent a conservation priority for the Park managers. There are nine locations distributed on siliceous soils both in the northwestern and southern slopes of the mountain range. In general, they are associated to major river valleys and within an altitudinal range of 1200 to 1900 m.a.s.l. (Table 1).

#### **2.2 Monitoring forest ecological status with EVI**

Our monitoring approach was based on the characterization of ecosystem functional attributes derived from the seasonal dynamics of the Enhanced Vegetation Index (EVI). The EVI calculates the normalized difference in reflectance between the red light that is absorbed in photosynthesis and the strong reflection of near infra-red light caused by the cell structure of the leaves. It also includes a third wavelength (blue) that is used to correct the influence of the atmosphere and the soil. EVI is defined according to equation 2 (Huete et al., 1997).

$$EVI = G \frac{NIR - R}{NIR + C\_1R - C\_2B + L} \tag{2}$$

Satellite-Based Monitoring of Ecosystem Functioning in Protected Areas:

1982). Significant trends were considered with p-values < 0.05.

*(Maximum activity)*

Maximum EVI (MAX)

0

2010) and G. Baldi from http://lechusa.unsl.edu.ar.

0.2

Relative

**2.3 Statistical analyses** 

*(Seasonality)*

Range or EVI\_sCV

0.4

0.6

0.8

1

**EVI**

Recent Trends in the Oak Forests (*Quercus pyrenaica* Willd.) of Sierra Nevada (Spain) 359

To explore the existence of inter-annual trends of ecosystem functioning during the 2001- 2009 period in the oak forests of Sierra Nevada, we followed the methodology suggested by Alcaraz-Segura et al. (2009b). In addition to the evaluation of long-term trends of the EVI\_mean, we also evaluated the existence of significant trends within each of the 23 images (16-day periods) of the year by means of the Mann-Kendall trend test, a non-parametric trend test robust against non-normality, heterocedasticity, outliers, and serial dependence. For each pixel, we obtained the slope of the trends through the Sen's Method (Hirsch et al.,

Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec

Date of maximum EVI (DMAX) *(Phenology)*

Fig. 2. Functional attributes of the EVI seasonal curve related to ecosystem primary production, seasonality, and phenology. EVI\_mean: EVI annual mean, EVI\_sCV: EVI seasonal Coefficient of Variation (SD/EVI\_mean), MAX: Maximum EVI annual value, MIN: Minimum EVI annual value, DMAX: Date in which is reached the maximum EVI value, DMIN: Date in which is reached the minimum EVI value. These attributes have a clear biological meaning, the EVI\_mean is an indicator of the fraction of the radiation used by plants and net primary productivity, MAX and MIN are indicators of the maximum and minimum photosynthetic activity, EVI\_sCV is one indicator of seasonality of carbon gains, and DMAX and DMIN are indicators of phenology. Image modified from (Cabello et al.,

To evaluate whether there exist differences in the EVI attributes among the nine oak woods studied in Sierra Nevada, we performed analysis of variance (ANOVA) only when either raw or transformed attributes fulfilled the necessary parametric requirements of normality and homoscedasticity. To reach normality, for EVI\_mean we applied a natural Logarithm (Ln) transformation (Shapiro-Wilk, W=0.990, p=0.266, n=177; Levene's Test F=0.474,

Date of minimum EVI (DMIN) *(Phenology)*

Minimum *(Minimum* EVI (MIN) *activity)*

Area under the EVI curve (EVI\_mean) *(Annual primary production)*

> **Time (months)**

Where NIR, R and B represent the reflectance in the near infrared, red, and blue wavelengths, C1 (6) and C2 (7.5) are coefficients of atmospheric resistance, G (2.5) is the gain factor, and L (1) is a soil correction factor.

Fig. 1. Distribution of the Pyrenean oak forest patches (*Quercus pyrenaica*) in the Sierra Nevada National Park (southeastern Spain). Forest patches are named according the river basin where they are located: Alhama, Genil, Monachil, Dílar, and Dúrcal, in the northern slope; and Chico, Soportújar, Poqueira, and Trevélez in the southern slope.

Our approach uses satellite images of the Enhanced Vegetation Index captured by the Moderate Resolution Imaging Spectroradiometer (MODIS) sensor onboard the Terra satellite from 2001 to 2009 (Product MOD13Q1). These images have a temporal resolution of 16 days (23 images per year) and a spatial resolution of 231x231 m. We used the Quality Assessment information to filter out low quality data, submitting images to a purification process which removes those pixels affected by high aerosol content, clouds, snow, shadows, and water. From this dataset, we first calculated the 9-year mean EVI seasonal curve for each oak forest site (Figure 1). For this, we only used pixels with more than 75% of their surface occupied by oak woods. Then, the following descriptive attributes of the ecosystem functioning were derived (Figure 2): The EVI annual mean (EVI\_mean), an estimator of primary production; the EVI seasonal (or intra-annual) coefficient of variation (EVI\_sCV), an indicator of seasonality of carbon gains; the EVI maximum (MAX) and minimum (MIN) values, indicators of the maximum and minimum photosynthetic capacities respectively; and the dates when the maximum (DMAX) and minimum (DMIN) EVI values are reached, two descriptors of the phenology of vegetation greenness. These attributes are widely used and have clear biological meanings (Pettorelli et al., 2005; Alcaraz-Segura et al., 2009a).

To explore the existence of inter-annual trends of ecosystem functioning during the 2001- 2009 period in the oak forests of Sierra Nevada, we followed the methodology suggested by Alcaraz-Segura et al. (2009b). In addition to the evaluation of long-term trends of the EVI\_mean, we also evaluated the existence of significant trends within each of the 23 images (16-day periods) of the year by means of the Mann-Kendall trend test, a non-parametric trend test robust against non-normality, heterocedasticity, outliers, and serial dependence. For each pixel, we obtained the slope of the trends through the Sen's Method (Hirsch et al., 1982). Significant trends were considered with p-values < 0.05.

Fig. 2. Functional attributes of the EVI seasonal curve related to ecosystem primary production, seasonality, and phenology. EVI\_mean: EVI annual mean, EVI\_sCV: EVI seasonal Coefficient of Variation (SD/EVI\_mean), MAX: Maximum EVI annual value, MIN: Minimum EVI annual value, DMAX: Date in which is reached the maximum EVI value, DMIN: Date in which is reached the minimum EVI value. These attributes have a clear biological meaning, the EVI\_mean is an indicator of the fraction of the radiation used by plants and net primary productivity, MAX and MIN are indicators of the maximum and minimum photosynthetic activity, EVI\_sCV is one indicator of seasonality of carbon gains, and DMAX and DMIN are indicators of phenology. Image modified from (Cabello et al., 2010) and G. Baldi from http://lechusa.unsl.edu.ar.

### **2.3 Statistical analyses**

358 International Perspectives on Global Environmental Change

Where NIR, R and B represent the reflectance in the near infrared, red, and blue wavelengths, C1 (6) and C2 (7.5) are coefficients of atmospheric resistance, G (2.5) is the gain

factor, and L (1) is a soil correction factor.

**Alhama**

**Genil**

**Monachil**

**Dílar**

**Dúrcal**

**Suportújar**

**Poqueira**

**Trevélez**

Fig. 1. Distribution of the Pyrenean oak forest patches (*Quercus pyrenaica*) in the Sierra Nevada National Park (southeastern Spain). Forest patches are named according the river basin where they are located: Alhama, Genil, Monachil, Dílar, and Dúrcal, in the northern

Our approach uses satellite images of the Enhanced Vegetation Index captured by the Moderate Resolution Imaging Spectroradiometer (MODIS) sensor onboard the Terra satellite from 2001 to 2009 (Product MOD13Q1). These images have a temporal resolution of 16 days (23 images per year) and a spatial resolution of 231x231 m. We used the Quality Assessment information to filter out low quality data, submitting images to a purification process which removes those pixels affected by high aerosol content, clouds, snow, shadows, and water. From this dataset, we first calculated the 9-year mean EVI seasonal curve for each oak forest site (Figure 1). For this, we only used pixels with more than 75% of their surface occupied by oak woods. Then, the following descriptive attributes of the ecosystem functioning were derived (Figure 2): The EVI annual mean (EVI\_mean), an estimator of primary production; the EVI seasonal (or intra-annual) coefficient of variation (EVI\_sCV), an indicator of seasonality of carbon gains; the EVI maximum (MAX) and minimum (MIN) values, indicators of the maximum and minimum photosynthetic capacities respectively; and the dates when the maximum (DMAX) and minimum (DMIN) EVI values are reached, two descriptors of the phenology of vegetation greenness. These attributes are widely used and have clear biological meanings (Pettorelli et al., 2005;

slope; and Chico, Soportújar, Poqueira, and Trevélez in the southern slope.

**Chico**

Alcaraz-Segura et al., 2009a).

To evaluate whether there exist differences in the EVI attributes among the nine oak woods studied in Sierra Nevada, we performed analysis of variance (ANOVA) only when either raw or transformed attributes fulfilled the necessary parametric requirements of normality and homoscedasticity. To reach normality, for EVI\_mean we applied a natural Logarithm (Ln) transformation (Shapiro-Wilk, W=0.990, p=0.266, n=177; Levene's Test F=0.474,

Satellite-Based Monitoring of Ecosystem Functioning in Protected Areas:

0 0.1 0.2 0.3 0.4 0.5

0 0.1 0.2 0.3 0.4 0.5

0 0.1 0.2 0.3 0.4 0.5

0 0.1 0.2 0.3 0.4 0.5

0 0.1 0.2 0.3 0.4 0.5

beginning and the end of the growing season.

**EVI**

1 Jan.

17 Jan.

2 Feb.

18 Feb.

6 Mar.

22 Mar.

7 Apr.

23 Apr.

9 May.

25 May.

10 Jun.

26 Jun.

Fig. 3. EVI seasonal dynamics (left Y axis, in gray) and 2001-2009 EVI trends (right Y axis, in black) for oak forests in the northern slope of Sierra Nevada. The horizontal "zero-trend" line shows the absence of significant trends. The two vertical dotted gray lines show the

12 Jul.

28 Jul.

13 Aug.

29 Aug.

14 Sep.

30 Sep.

16 Oct.

1 Nov.

17 Nov.

3 Dec.

19 Dec.

**EVI**

1 Jan.

17 Jan.

2 Feb.

18 Feb.

6 Mar.

22 Mar.

7 Apr.

23 Apr.

9 May.

25 May.

10 Jun.

26 Jun.

**e) Dúrcal**

12 Jul.

28 Jul.

13 Aug.

29 Aug.

14 Sep.

30 Sep.

16 Oct.

1 Nov.

17 Nov.

3 Dec.

19 Dec.

**EVI**

1 Jan.

17 Jan.

2 Feb.

18 Feb.

6 Mar.

22 Mar.

7 Apr.

23 Apr.

9 May.

25 May.

10 Jun.

26 Jun.

**d) Dílar**

12 Jul.

28 Jul.

13 Aug.

29 Aug.

14 Sep.

30 Sep.

16 Oct.

1 Nov.

17 Nov.

3 Dec.

19 Dec.

**EVI**

1 Jan.

17 Jan.

2 Feb.

18 Feb.

6 Mar.

22 Mar.

7 Apr.

23 Apr.

9 May.

25 May.

10 Jun.

26 Jun.

**c) Monachil**

12 Jul.

28 Jul.

13 Aug.

29 Aug.

14 Sep.

30 Sep.

16 Oct.

1 Nov.

17 Nov.

3 Dec.

19 Dec.

**EVI**

1 Jan.

17 Jan.

2 Feb.

18 Feb.

6 Mar.

22 Mar.

7 Apr.

23 Apr.

9 May.

25 May.

10 Jun.

26 Jun.

**b) Genil** 

12 Jul.

28 Jul.

13 Aug.

29 Aug.

14 Sep.

30 Sep.

16 Oct.

1 Nov.

17 Nov.

3 Dec.

19 Dec.

Recent Trends in the Oak Forests (*Quercus pyrenaica* Willd.) of Sierra Nevada (Spain) 361

**a) Alhama** 






**EVI Trends**

**EVI Trends**

**EVI Trends**

**EVI Trends**

**EVI Trends**

p=0.873, n=177) and for EVI\_sCV a Box-Cox transformation (Shapiro-Wilk, W=0.983, p=0.031, n=177; Levene's Test F=1.951, p=0.055, n=177). The slight but not significant deviation from normality for the EVI\_sCV data did not affect results. For those attributes that even transformed did not fulfill normality (MAX, MIN, DMAX, and DMIN), the analysis was conducted using the non-parametric Kruskal-Wallis test. To determine which groups significantly differed from each other, we used multiple *post hoc* comparisons, using the Tukey test for EVI\_mean and EVI\_sCV, and the Bonferroni test for MAX, MIN, DMAX, and DMIN. See Figure 5.
