Selective logging changes forest phenology in the Brazilian Amazon: Evidence fromMODIS image time series analysis
Alexander Koltunov a,⁎, Susan L. Ustin a, Gregory P. Asner b, Inez Fung c
a Center for Spatial Technologies and Remote Sensing, University of California, Davis, The Barn, One Shields Avenue, Davis, CA 95616, USAb Department of Global Ecology, Carnegie Institution, Stanford, CA 94305, USAc Department of Earth and Planetary Science, University of California, Berkeley, CA 94720, USA
We present a large-scale study of the relationships between selective logging and forest phenology in theBrazilian Amazon. Time-series analysis of MODIS satellite data of selectively logged forests in Mato Grosso,Brazil, shows that relatively low levels (5–10%) of canopy damage cause significant and long-lasting (morethan 3 years) changes in forest phenology. Partial clearing slows forest green-up in the dry season,progressively dries the canopy, and induces overall seasonal deficits in canopy moisture and greenness. Givenlarge and increasing geographic extent of selective logging throughout Amazonia, this phenologicaldisturbance has a potential to impact carbon and water fluxes, nutrient dynamics, and other functionalprocesses in these forests.
Selective logging or partial forest removal in the Brazilian Amazonwas recently shown in analyses of Landsat Enhanced ThematicMapper Plus (ETM+) data at high spatial resolution (30-m pixels)to be occurring at rates of about 12,000–20,000 km2 per year (Asneret al., 2005), thus indicating the central role of selective logging intropical forest disturbance. Selective timber harvests in the State ofMato Grosso are reported to range from 10,000 to 40,000 m3 per km2
(Monteiro et al., 2004; Souza et al., 2005b), following clusteredremoval of trees, construction of skid trails, landings, and roads thatfrequently damage remaining forest (Matricardi et al., 2005; White,1994; Asner et al., 2004). Although selective harvests have far lessimpact on forest processes than deforestation, selectively logged sitesexperience higher rates of forest fires (Cochrane, 2003; Matricardiet al., 2005), tree fall (Young and Hubbell, 1991), changes inmicroclimate (Ray et al., 2005), soil compaction and erosion (Olanderet al., 2005), among other impacts.
Removal of merchantable species during selective logging changesforest composition and structure in a way that may alter forestfunction and biosphere–atmosphere exchanges of carbon dioxide,water vapor, and energy. Any such changes in forest function might beexpressed by changes in forest phenology — the seasonal cycle ofcanopy development and physiological processes. Climate seasonalityin parts of the Brazilian Amazon is well pronounced, with a dry season
).
l rights reserved.
often lasting 4–5 months and demonstrated linkages to forest foliageand greenness (Huete et al., 2006; Xaio et al., 2006), and to forestfunctional properties (Myneni et al., 2007).
1.1. Related work and problem statement
Deforestation is the conversion of forest to another land-covertype, often with profound impacts on phenology of the remaining ornew (e.g. pasture) vegetation (Pongratz et al., 2006). Little informa-tion exists onwhether land uses that leave the forest standing, such asselective logging, might have an impact on phenological dynamics.Few studies have analyzed the effects of anthropogenic disturbance onforest phenology. Herrerias-Diego et al. (2006) found that fragmenta-tion shifted the flowering period of a perennial tree species by 15–20 days for at least three years after the disturbance. Such a 2–3 weekdelay could impact fruit maturation.
Previous multitemporal analyses of selectively logged areas havebeen limited to inter-annual differences (Stone and Lefebvre, 2000;Souza and Baretto, 2000; Souza et al., 2003; Souza et al., 2005a; Souzaet al., 2005b; Souza and Roberts, 2005) or to field measurementscollected at a few plots (Schulze and Zweede, 2006). These studieshave been primarily focused on improving selective logging detectionalgorithms, in which detection of canopy gaps and analysis of theirdynamics play a critical role. Several satellite-based studies indicatethat foliage cover returns in just 1–3 years following logging (Stoneand Lefebvre, 2000; Souza and Baretto, 2000; Souza et al., 2003; Souzaand Roberts, 2005; Souza et al., 2005b; Matricardi et al., 2005), mostlybased on inter-annual similarity of vegetation indexes (VIs) such as
2432 A. Koltunov et al. / Remote Sensing of Environment 113 (2009) 2431–2440
the Normalized Difference Vegetation Index (NDVI, (Tucker, 1979)).Souza et al., 2005b reported an increase in satellite-based fractionalcanopy greenness, primarily in the first year after logging, althoughsome effects persist longer for various types of logging practices(Matricardi et al., 2005). Others report recovery periods lastingdecades for regeneration of primary forest structure and leaf area(Feldpausch et al., 2005; Okuda et al., 2003). Asner et al., 2004recently provided a quantitative assessment of the gap closure ratesover large areas in the Brazilian Amazon and showed that approxi-mately 50% of the canopy openings caused by logging are closedwithin one year of regrowth following timber harvests. This recoveryapplies only to foliar cover, and not to wood biomass, which doesindeed take many decades to recover (Keller et al., 2004).
Rapid closure of relatively small canopy gaps following selectivelogging may appear similar to natural regeneration phenomena inthese forests. As a result, the phenological trajectories related to forestfunction and biosphere–atmosphere exchange may seem unaffected.In turn, the absence of significant phenological changes wouldsupport the claim that selective logging is effective at maintainingbackground conditions of tropical forests.
As a counterargument, it is not immediately obvious, if at all true,that selective logging alters biospheric processes in a manner limitedto gap closure processes. Indeed, even with gaps completely closed,changing forest composition and age class structure may impact forestphenology either alone or in combination with other environmentalfactors (Young and Hubbell, 1991; Ray et al., 2005). Soil compactionand erosion could also affect water and nutrient availability that affectphenological processes independent of crown closure. In other words,analysis of gap dynamics and analysis of phenological changes in thelogged forests are generally not equivalent.
Given the prevalence of selective logging in Amazonian forests,and the uncertainty about forest function following timber harvest,we sought to answer the basic question:
1. Does selective logging change the phenological dynamics of theforest in the Brazilian Amazon?
And, if it does:
2. What are the characteristic features of these changes?3. How long do the differences remain significant?4. Are the changes pertinent only to highly damaged forests and
insignificant if the initial canopy damage was low?
To date, none of these questions have been addressed.
1.2. Phenological change analysis
Small-scale studies, experimental or observational, may provideuseful clues about possible phenological consequences of selectivelogging, and are necessary to develop testable (alternative) hypoth-eses. However, conclusions from small-scale studies cannot be readilygeneralized to regional or continental scales, because variability oflocalized environmental factors is underrepresented in small areas.Furthermore, tropical forest phenology is a nearly periodic dynamicprocess with seasonal variability much larger than inter-annualvariability (Huete et al., 2006). Therefore, analysis of a singleobservation of forest greenness per year or per season cannotcharacterize phenology and will not lead to reliable inference aboutthe magnitude, direction, or even the presence of phenologicalchanges.
To detect and characterize the phenological changes associatedwith selective logging requires frequent observations over largegeographic areas, necessitating the use of low spatial resolution datafrom weather satellites. However, at the low spatial resolution ofweather satellite data, the spectral signal of small-scale forestclearings is subtle (Fig. 1), making identification of these impactsand separation from background forests extremely difficult. For
example, using the MODerate Resolution Imaging Spectroradiometer(MODIS) satellite data, Morton et al., 2005 were able to detectdeforestation polygons that were greater than three 250-m pixels(~19 ha) in Mato Grosso, Brazil, but were unsuccessful in detectingclear-cuts of less than half-pixel size. Yet, in MODIS imagery, selectivelogging events appear at far smaller scales, and in a more diffusepattern than clear-cutting, making them even less distinct from intactforests.
To quantify damages induced by selective logging one also needsan index representing logging intensity. Addressing this issue, Asneret al. (2006) developed a measure of canopy damage due to selectivelogging by converting remotely sensed subpixel fractional changes incanopy cover (derived from Landsat 7 data) into spatial estimates offorest canopy gap fraction. Because forest gaps are larger and morenumerous in areas disturbed by logging (i.e. felling gaps, skid trails,and logging decks), this measure, termed Area-Integrated GapFraction (AIGF), provides a consistent estimate of canopy damagewithin logged forests. Higher AIGF values indicate amore open canopyand thus greater disturbance throughout the forest. More details onthe AIGF measure are given below in Section 2.2.
Finally, to study the effect of selective logging on phenologicalbehavior of the forest, one needs to know which forests wereselectively logged or clear-cut in the years following the initialselective logging events. These forests should be removed fromsubsequent analysis to avoid confounding the results. Currently,selective logging data over the large area under study are available inthe form of AIGF for year 2000 and in the form of binary (0/1) maps ofselective logging events for years 2001 and 2002, thereby limiting thefeasible temporal range of phenology analysis to these three years.
1.3. Study objectives and outline
The objectives of this study are to uncover and quantify theregional-scale effects of selective logging on forest phenology. Westudied phenological consequences of logging that occurred between1999 and 2000 in the State of Mato Grosso, Brazil (Fig. 1), using a timeseries of MODIS multispectral 1-km resolution imagery combinedwith AIGF estimates derived from Landsat 7. Thus, any phenologicalimpacts detected at the individual pixel scale represent changes inphenological responses that are integrated over 1-km2 area of loggedforest, including not only the actual tree canopies, but also changes insoil, understory, and other background effects. These phenologicalimpacts of logging were computed and analyzed for the same timeinterval each year: May 8–August 13 in 2000, 2001, and 2002,representing the dry season in Mato Grosso.
2. Materials and methods
2.1. Study area and satellite image sequence
Our study area is located in the northern part of the State of MatoGrosso, and confined within a bounding rectangle with cornercoordinates: −8.8 W, −61.53 S (upper-left); −13.75 W, −50.34 S(lower-right). The area (Fig. 1) of nearly 670,000 km2 is characterizedby a large number of selective logging events that occurred between1999 and 2000, with total logged area exceeding 13,000 km2 (Asneret al., 2005), whichwas considered a high logging year in terms of areaof timber harvest, thus providing a large number of pixels forsubsequent phenological analysis.
MODerate Resolution Imaging Spectroradiometer (MODIS) is a keyinstrument onboard the Terra and Aqua polar-orbiting satellites. EachMODIS instrument provides near-daily coverage of the entire Earth'ssurface, acquiring data in 36 bands in visible and infrared spectraldomain under different solar-view geometries at the spatial resolutionof 250 to 1000 m, depending on the wavelength. The bi-directionalreflectance distribution function (BRDF) retrieval method used in the
Fig. 1. False-color composite image and location of the study scene whose area is nearly 670,000 km2. The MODIS band wavelengths [µm] used for composition: red: 2.1, green: 8.6,blue: 6.5. Vegetation appears green, water — black or dark blue, non-forested areas — pink, deforestation — greenish gray to light pink. The zoom window (upper-left corner)illustrates the apparent spectral similarity of selectively logged and unlogged forests at 1-km scale.
2433A. Koltunov et al. / Remote Sensing of Environment 113 (2009) 2431–2440
MODIS BRDF products (Schaaf et al., 2002) employs a constrainedsemi-empirical phenomenological BRDF model (Roujean et al., 1992)whose time-invariant parameters are computed for each pixel. Thismethod uses multi-date multi-angular views acquired under differentillumination geometries over a 16-day period. In this way, theresulting NBAR (Nadir BRDF Adjusted Reflectance) products (Strahleret al., 1999), although delivered each 16 days, integrate the reflectanceinformation from all cloud-free observations that were acquired on anear-daily basis, over the 16-day period. The NBAR images areadjusted for aerosol, atmospheric, and bi-directional effects (Strahleret al., 1999; Schaaf et al., 2002) and masked for clouds and cloudshadows (Strahler et al., 1999). We compiled 65 MODIS NBAR imagesfrom Terra spanning the interval fromMarch 15, 2000 to December 31,2002, at 16-day time steps and at 1-km ground resolution.
2.2. Area-integrated gap fraction (AIGF)
2.2.1. Estimating forest canopy gap fraction in Landsat imageryForest canopy gap fraction, or the fractional canopy cover in the
upward-pointing hemisphere from any given ground location, haslong served as a measure of canopy structure in forest ecosystems.Canopy gap fraction largely determines photosynthetic rates, canopyenergy and water balance, primary production, mammal and insectdynamics, and even the probability of fire. It is therefore highlyadvantageous to convert a remotely sensed radiometric measurementof fractional photosynthetic vegetation (PV) cover to the traditionalcanopy gap fraction measure. The two measures are not the same.First, PV fraction is a planar metric, whereas canopy gap fraction ishemispherical in nature. Moreover, strong adjacency effects between
satellite pixels (caused by interpixel light scattering) result in anonlinear component to vegetation mixture modeling with multi-spectral data. This effect is maximum at forest gap values greater than≈85% (Asner et al., 2005).
Asner et al. (2005) developed a set of equations relating PV coverfraction derived from spectral mixture analysis of Landsat imagery tofield-based measurements of forest canopy gap fraction. Asner et al.(2006) improved the PV–gap fraction relationship, this time via amore extensive comparison of PV fractions from spectral mixtureanalysis to field-based forest canopy gap fractions collected across awide range of low-, medium- and high-impact logging sites. Themethod provides the means to convert Landsat reflectance imagery tofractional PV cover to forest canopy gap fraction. An example for onearea within the State of Mato Grosso, Brazil, is shown in Fig. 2.
2.2.2. AIGFFollowing derivation of canopy gap fraction maps from Landsat
imagery, a 1×1-km window can be passed over the gap maps tocalculate the area-integrated gap fraction (AIGF):
AIGF =P
fin
; ð1Þ
where fi is each forested pixel's gap fraction value within each 1-km2
area, and n is the total number of pixels within that area. Becauseforest gap is greater in areas having been disturbed by logging (i.e.,felling gaps, skid trails, and logging decks), the AIGF provides aconsistent estimate of canopy damage within logged forest, where a
Fig. 2. Forest canopy gap fraction derived from spectral mixture analysis of Landsatimagery (from Asner et al., 2005).
2434 A. Koltunov et al. / Remote Sensing of Environment 113 (2009) 2431–2440
higher AIGF value indicates a more open canopy and thus moredisturbance throughout the forest as a whole.
Asner et al. (2006) used maps of AIGF derived from Landsat toquantify post-logging forest canopy closure. The means and standarddeviations of AIGF used in that analysis were calculated from over180,000 pixels per AIGF class. Mapping AIGF is a straightforwardapproach to quantifying the recovery of logged forest over time.
2.3. Estimation and modeling of phenological impacts
Let AIGF0 denote AIGF of selective logging that occurred between1999 and 2000, which was observed in 2000 and mapped at 1-kmscale by Asner et al., 2006. We use AIGF0 as a logging intensitymeasure in the analyses described below.
The INPE PRODES (Program for Monitoring Deforestation in theBrazilian Amazon) deforestation map for year 2002, the 1-km scaleAIFG0 map, and two 30-m resolution maps of selective logging eventsoccurred between 2000–2001 and 2001–2002, respectively (Asneret al., 2005), allowed us to extract 74,432 pixels suitable for ouranalysis. These were pixels that were not clear-cut by 2003 (i.e. theywere labeled “forest” in the PRODESmap for year 2002), nor selectivelylogged (AIGF b0.05% ) between year 2000 and 2003, based onmaps ofselective logging events (Asner et al., 2005). To facilitate the analysis,these pixels were split into nine disjoint groups, G0, G1,…, G8, withrespect to AIGF0. Thus, for example, group G8 is a logged forest withAIGF0 70 to 80%, and group G0 represents forests with no knownhistory of logging (gap fraction b0.05% ) during 1999 to 2003. The firstfour columns of Table 1 provide detailed information about the groups.For convenience, we refer to group G0 as unlogged or intact forest.
In each pixel, we used MODIS NBAR images to calculate the timeseries of two vegetation indexes, the Enhanced Vegetation Index [EVI,(Huete et al., 1994)] and the Normalized Difference Water Index[NDWI, (Gao, 1996)], which estimate vegetation greenness andmoisture content, respectively. These indexes use MODIS reflectancebands ρ in red (band 1), blue (band 3), near infrared (NIR, band 2),and shortwave infrared (SWIR, band 5), according to the definitionsgiven below:
EVI = G ×ρNIR − ρRED
ρNIR + C1ρRED − C2ρBLUE + LNDWI =
ρNIR − ρSWIR
ρNIR − ρSWIRð2Þ
where, G=2.5, C1=6.0, C2=7.5, and L=1.0. The EVI reduces theeffect of aerosol influences in the red band and has a higher dynamicrange than the more traditional NDVI over the more humid forestedareas (Huete et al., 2002). EVI has also been shown to be well suitedfor monitoring phenological trajectories (Xaio et al., 2005) of tropicalforests.
2.3.1. Impact trajectory estimationSatellite-based estimation of the phenological impact of a
disturbance event entails estimation of the spectral signal thatwould have been observed at a spatial location s at time t, if therewas no disturbance. In general, two approaches can be followed:temporal and contextual. The temporal approach requires spectralmeasurements before the disturbance events. Due to unavailability ofMODIS data before 2000, the only approach to follow was thecontextual one. To estimate the intact forest signal at location s andtime t, we used neighboring pixel data at the same time t in a 20-kmwindow centered at s. For each vegetation index, the phenologicalimpact time series, ΔV (s, t), was computed as the difference betweenthe observed VI in the window center and the mean VI of the validneighbors, i.e.:
ΔV s; tð Þ = V s; tð Þ− μ s; tð Þ: ð3Þ
Valid neighbors are non-clouded pixel-members of G0 of the samevegetation type, according to the vegetation map (IBGE, 2005), as thecentral pixel.
The relative error of estimation, E=E(s, t), is determined by theheterogeneity of the spectral response within a window:
E sð Þ = V sð Þ− μ sð ÞV sð Þ =
ΔV sð ÞV sð Þ ; for each t: ð4Þ
The errors for unlogged pixels are known and can be used to assessthe estimation accuracy over selectively logged pixels. Fig. 3 displaysthe time series of the relative RMSE of EVI and NDWI estimation overunlogged forests for years 2000 to 2002, and demonstrates that thecontextual analysis over valid neighbors is a reasonably accurateapproach, with temporal variability of error from 3 to 5% in the dryseason. The reliability of the error (4) in impact estimation dependson the number of valid pixels inside the window, among other factors.As selectively logged events exhibit obvious spatial clustering and areoften located next to deforested areas, choosing a small window sizewould bias the estimated signal, µ(s), of unlogged forest. In turn, therelatively large window size (20 km) increases the number of validpixel-neighbors per pixel in question. In general, suboptimal windowsizes tend to increase the error E(s) of impact estimation and/or thevariance of this error. Future work may improve the estimation of theintact forest signal by analyzing the local autocorrelation field perpixel.
Considering impacts defined by (3) rather than the VIs themselves,provides a means to minimize confounding effects of environmentalfactors correlating with geographic location.
2.3.2. Parametric seasonal modeling of greenness and moisture impactsFor a given VI and a given time interval [a, b] representing a season,
we introduce a seasonal phenological model of the impact time series(defined by Eq. (3)) of a pixel:
where the basis functions e0(t),…,eK(t) are discrete orthogonalpolynomials on [a, b] (as specified in Appendix I). The coefficients ofthis model, α0,…, αK, characterize the mean amplitude (α0) and thedynamic pattern (α1,…, αK) of the changes in the forest VI over the
Table 1Gap fraction groups and upper bounds for p-values (u.b.p.v.) of significance tests.
Shownoutside theparentheses areu.b.p.v. of the permutation test for significanceof negative difference in themean impact parameter ofGk and that ofG0, (k≠0). Shown inparentheses areu.b.p.v.of the t-test that themean impact parameter of Gk is zero, against the alternative hypothesis that themean impact parameter of Gk is negative (k≠0). For G0, u.b.p.v. represent the two-sided one-sample t-tests. Upper bounds not exceeding 10−7 are displayed as zeros.
2435A. Koltunov et al. / Remote Sensing of Environment 113 (2009) 2431–2440
season. Error of representation (5) is assumed to have a normaldistribution. Because the polynomials are orthogonal, for each pixel,the least squares estimators of the coefficients are uncorrelated andindependent of the approximation order K. Consequently, theyrepresent independent features of the phenological impact dynamics.
We consider the first two coefficients of this model, α0 and α1, andterm them as follows:
• when using MODIS EVI:- α0: the Effective Intra-Seasonal Greenness Impact (SGI).- α1: the Effective Intra-Seasonal Greenness Impact Trend (SGIT orthe SGI trend)
• when using MODIS NDWI:- α0: the Effective Intra-Seasonal Moisture Impact (SMI).
Fig. 3. Times series of relative RMSE for contextual estimation of MODIS EVI and NDWIof unlogged forests. Time step — 16 days.
Thus, SGI and SMI represent overall (mean) impacts on [a, b], whileSGIT and SMIT characterize the linear trend of the impact on [a, b].Altogether, we call these metrics “impact parameters”. Eq. (3) impliesthat for unlogged forests the expected value of each impact parameteris zero.
In general, the choice of the time interval [a, b] influences theshape of the phenological trajectories and therefore the parameters ofthe impact model (5). The multitemporal sequence of satellite imagesused in this study spanned three dry seasons for years 2000, 2001, and2002 in Mato Grosso. In each year, we selected the same interval, May8–August 13, with seven NBAR images available. This intervalrepresents a median dry season in Mato Grosso (Table 2) and ischaracterized in the 2000–2002 period by the least cloudy MODISobservations during the year. If a longer interval is considered, the
Table 2Dry season timing, according to monthly averaged temperature and precipitation datain Mato Grosso.
Weatherstation
Geographiccoordinates
Begin (Julian Day/30) End (Julian Day/30)
Cuiaba 15°39′ S, 56°06′ W 4.6 8.6Diamantino 18°14′ S, 43°36′ W 4.6 8.0Alto Tapajos 7°21′ S, 57°31′ W 5.2 6.7Sao Felix Xingu 6°38′ S, 51°58′ W 5.4 7.3Conceicao DoAraguaia
8°15′ S, 49°12′ W 4.7 8.3
Sangradonso 15°38′ S, 53°54′ W 3.4 8.6Aragarcas 15°54′ S, 52°14′ W 4.0 8.2Porto Nacional 10°31′ S, 48°43′ W 4.4 8.5Median 4.60 8.25
The source data can be obtained at: http://www.globalbioclimatics.org/.
2436 A. Koltunov et al. / Remote Sensing of Environment 113 (2009) 2431–2440
intra-seasonal impact analysis becomes less reliable, because after theresumption of the rains, the phenological behavior of the vegetationbecomes more heterogeneous and therefore the spectral featurescharacterizing the dry season become mixed with spectral features ofthe wet season. The wet season was not followed due to loss of datafrom frequent cloud cover.
Foreachyear: 2000,2001, and2002,wemodeledandtestedonly thosepixels which were not clouded or had missing values in any of the NBARimages of the seasonal interval, and which also had at least 25 validneighbors in each NBAR image. The intra-seasonal impact parameters forthe three dry seasons were computed by fitting themodel (5) in the leastsquares sense using the orthogonality property of the basis functions.
2.4. Testing significance of phenological impact
For each impact parameter, the group-conditional distributionshave similar but not identical shape and standard deviations (with std[G0]/std[Gk] ranging from 0.96 to 1.11) and they do not exactly follow anormal distribution, which may affect accuracy of the standard t-testp-values. In these circumstances, in addition to the left-sided one-sample t-tests, we applied non-parametric permutation significancetests (Good, 2005) of the differences in means and in 10%-trimmedmeans for impact parameters of each selective logging class (Gk) andfor the unlogged forest (G0). The alternative hypothesis was that thedifferences were negative. The same number of permutations, 104,were used in all two-sample tests. The upper bounds for the one-sample t-test p-values are shown in parentheses in Table 1.
Although providing large samples, the use of unlogged neighbor-ing pixels relies on the implicit assumption that the probability of afuture logging event at a forest pixel does not depend on the VI valuesof that pixel. To assess the influence of this assumption on the impactparameter analysis, we compared the parameters of the forest loggedin the years following the year of comparison (6294 pixels — in 2001and 5372 pixels — in 2002) with the unlogged forest (group G0). Atwo-sample permutation significance test with 104 replicationsrevealed the significance of the negative differences between theimpact trends of the forest to be logged in subsequent years and thoseof the unlogged forest (Table 3). This indicated the presence ofconfounding variables that precluded establishing a cause–effectrelationship between selective logging and the observed phenologicaldifferences. Therefore, to enable causal inferences the class-conditionaldistributions of each impact parameter for classes G1,…,G8 werecorrected for the corresponding bias before conducting the testssummarized in Table 1. The bias of the impact parameter αk (k=1,2)of the dry season of year Y (Y=2000, 2001) was computed as thedifference µ0,L(Y)−µ0(Y), where µ0,L(Y) is the mean value of αk for thepixels that were unlogged in year Y but were logged in year (Y+1); andµ0(Y) is the mean value of αk for the pixels with no known historyof logging (i.e. from group G0). The parameter bias for 2002 wasestimated by averaging the biases observed in 2000 and 2001.
Table 3Results of the permutation tests for significance of selective logging occurrence in areaswith negative SGI trend and SMI trend.
Groups Pixelstested
Alternative hypothesis Upper bound forp-value in 2000
Upper bound forp-value in 2001
Logged after2000 only
6294 Mean[SGIT]bmean[SGIT(G0)]
b10−4 (0) Not applicable
Mean[SMIT]bmean[SMIT(G0)
b10−4 (0)
Logged after2001 only
5372 Mean[SGIT]bmean[SGIT(G0)]
b10−4 (0) b10−4 (0)
Mean[SMIT]bmean[SMIT(G0)
b0.001 (0) b10−4 (0)
Shown in parentheses are the upper bounds for p-values of the t-test that the meanimpact trend of the corresponding group is zero, against the alternative hypothesis thatit is negative. Upper bounds not exceeding 10−7 are displayed as zeros.
The upper bounds for all two-sample test p-values in all tableshold also for 10%-trimmed means. They are shown outside theparentheses in the tables.
3. Results and discussion
3.1. Phenological changes
The effects of selective logging on greenness and moisture areobserved, relative to the intact forest (Fig. 4), by the downwarddisplacement of the mean impact trajectories of selectively loggedforest classes. Also, the magnitude of differences between selectivelogging and intact forest markedly increases with time from thebeginning to the end of the dry season. All selectively logged classesexhibit a similar dynamic pattern of impacts over the dry season, inmarked contrast with the unlogged forest. As shown in studies ofother ecosystems (Cheng et al., 2008), EVI and NDWI follow similarbut not identical temporal trajectories, indicating that these indexesprovide some independent information about canopy functionality.
Our statistical analysis confirms that, in the dry season of each year,SGI, SMI, SGIT, and SMIT of selectively logged forests are significantlynegative and less than those of unlogged forests (Table 1), mostlyexhibiting a monotonic increase in magnitude as AIGF0 grows (Figs. 5and 6). The greenness impact parameters, SGI and SGIT, are significanteven to the 5–10% disturbance class level, and through the third yearafter selective logging.
With years following logging events, the magnitude of SGI and SMIdecreases at the observed rate of 52% (±11%) a year (Figs. 5 and 6A, B),consistent with field observations of gap closure rates by Asner et al.,2004. After three years, the overall VIs are nearly recovered to thevalues of neighboring unlogged forest, consistent with previousremote sensing observations (Stone and Lefebvre, 2000; Souza andBaretto, 2000; Souza et al., 2003; Souza and Roberts, 2005; Souza et al.,2005b; Matricardi et al., 2005). The continued significant negativedifferences in SGI and SMI three years after selective logging likelyindicate incomplete gap closure, altered canopy structure, speciescomposition differences, or increased vulnerability and higher post-logging mortality of the remaining trees (Young and Hubbell, 1991;Schulze and Zweede, 2006).
The selectively logged sites tended to occur in forest areas thatalready had negative and smaller SGIT and SMIT than the forests thatremained unlogged in 2002 (see Section 2.4). This difference waslikely caused by proximity to a road or sites logged in previous years,which were located within the 1-km2 area of the pixel. It has beenpreviously observed that these factors increase the probability oflogging in nearby locations. Nonetheless, after correction for biasesdue to pre-existing negative SGIT and SMIT, both parameters remainsignificantly negative for selective logging (Table 1). That is, the dataprovide strong evidence that selective logging changes the phenolo-gical pattern of forests by inducing relative drying and impeding thedry season green-up, a phenomenon recently reported for unloggedforests by Huete et al., 2006, Xaio et al., 2006, and Myneni et al., 2007.
Furthermore, not only do the SGI trend and SMI trend not recoverin the three year time interval, consistent with the gap closure rates,but additionally, there is no evidence that the impact trends areconverging to those of the intact forest. This can be seen in Fig. 6,which illustrates the directions of inter-annual change in the impactparameters. The lack of convergence of SGIT and SMIT to zero withyears (Fig. 6C, D) indicates that the combined aboveground andbelowground changes in forest structure are progressive or unrecov-ered on short multi-year time scales, especially for damaged forestswith more than 10% canopy gap in 2000.
Remarkably, the changes in phenological patterns of greenness aresignificant down to the lowest damage group (AIGF0 5–10%) compris-ing 54% of all year 2000 logging observations. This finding supports thehypothesis that spreading logging operations over large areas does not
Fig. 4. Temporal trajectories of intra-group means (solid lines) of the VI differences from the neighboring intact forest. Dashed lines depict mean±standard deviation. Wet seasondata are omitted.
2437A. Koltunov et al. / Remote Sensing of Environment 113 (2009) 2431–2440
entirely mitigate functional changes to forests following timberextraction. Given no indication that the impact trends approach zerowithin three years after logging, continued selective harvest practicesin this forest will accumulate low per-hectare damage with time.
3.2. Possible direct causes of the observed changes
Obviously, selective logging changes the observed forest phenolo-gical signal indirectly. We believe that the list of more direct causes ofthe observed changes varies in space and time, including but wouldnot be limited to: forest herbaceous understory and younger treesfilling the gaps, soil compaction and erosion, as well as changes inmicroclimate and nutrient cycling, each of which could produceenhanced drying of the selectively logged forest relative to its normalconditions. These phenomena have been previously reported to followforest clearing (Nobre et al., 1991; Ter Steege et al., 1995; McNabbet al., 1997; Asner et al., 2004; Ray et al., 2005), thereby providingindirect field-based support for mechanisms that could account forthe observed phenological changes. Perhaps the major factor is thereduction in the average root depth due to selective logging.
Herbaceous species and younger trees have limited or no access todeep subsurface water, which become increasingly important over thecourse of the dry season. Consequently, the ability of these forests tomaintain and increase (Huete et al., 2006) photosynthetic activity andevapotranspiration during the season may be progressively reduced,resulting in the observed changes (Fig. 4) in satellite-derived area-integrated phenological trajectories. Future field-based studies areneeded to complete the list of potential biological and physiologicalmechanisms and accurately estimate their relative contribution to theobserved impacts. Whatever the actual direct mechanisms are, theysignificantly alter forest phenology only after selective logging.
3.3. Insignificant factors
In addition, our study suggests which factors following selectivelogging are unlikely to individually explain our satellite observations.For instance, post-logging clear-cut activities (Asner et al., 2006), evenif missed by the PRODES deforestation maps, are unlikely to affectgreenness proportionally to AIGF0 (Figs. 5 and 6). Similarly, pheno-logical changes were significant with respect to both measures: group
Fig. 5. Sample scatter plots of impact parameter distributions of AIGF groups. Solid lines show within-group means. The considerable overlap of the within-class parameterdistributions illustrates a low impact-to-noise ratio, and therefore the difficulty of confirming the existence of the impact. Note that the vertical lines (mean overall intra-seasonalimpacts) shrink to zero by ~50% each year, following the canopy closure processes, whereas horizontal lines (mean trends of the impact) do not exhibit a similar inter-annualbehavior.
2438 A. Koltunov et al. / Remote Sensing of Environment 113 (2009) 2431–2440
mean and 10%-trimmed mean, and therefore may not be attributed tofires following selective logging or to any other low-probabilitycompounding disturbance events (e.g., defoliation by insects). Theeffect of varying subpixel shadow fraction in the satellite-derived
Fig. 6. Inter-annual dynamics of the mean impact parameters of selective logging classes. Noteach year, according to the canopy closure rates, whereas the magnitude of the mean impa
phenological changes is marginal, because it is corrected for in theMODIS NBAR images (Schaaf et al., 2002). Furthermore, canopy gapclosure is a monotonic process particularly occurring throughout thedry season, owing to increasing leaf area in the forest (Xaio et al.,
e that the magnitude of the mean overall intra-seasonal impacts (A, B) reduces by ~50%ct trends (C, D) does not exhibit a similar inter-annual pattern of behavior.
2439A. Koltunov et al. / Remote Sensing of Environment 113 (2009) 2431–2440
2005; Myneni et al., 2007), and thus contributing to a positive intra-seasonal temporal trend (Huete et al., 2006) in areal measurements offorest greenness, such as EVI. Therefore, had gap dynamics been theonly driver of post-logging phenological differences, a significantlynegative SGIT would not have been observed.
3.4. Spatial scale of inference
The existence of phenological impacts due to selective logging isconfirmed with high confidence (Table 1) over a territory comparablein size to the Brazilian Amazon.
For a given individual pixel area, selective logging acts as a stressagent, which (under absence of confounding factors) increases theprobability that this area of forest will undergo phenological changesin the directions described in Section 3.1. As a consequence of theindirect and probabilistic nature of this cause–effect relationship, wesee that for individual pixels the impact may be higher, lower, orabsent, depending on forest sub-groups, geographical or biophysicalvariables, different climate conditions, or other circumstances. Thesefactors influencing the observed VIs, the limitations of the contextualestimation of the impact trajectories, and as mentioned in Section 1.2,the low spatial resolution of the satellite data, contribute to the lowsignal-to-noise ratio (SNR) of the obtained impact data, especially forthe AIGF classes below 30% (Fig. 4). It is well known that a small SNRdoes not definitively imply the absence of the signal. What it doesimply is that the signal could be absent or it is poorly visible. That iswhy the large number of randomly selected SL events and therobustness of the test procedures have played a critical role inconfirming with high confidence that selective logging does changeforest phenology in the Brazilian Amazon.
Similar to VIs and AIGF, the impact parameters used in this studyare effective quantities, meaning that they integrate the impact overthe pixel area. Therefore, quite obviously, the phenological impactestablished by this paper is valid for forests, i.e. areas with a highdensity of trees. Whether or not the remaining trees themselves tendto dry out due to selective logging or only understory species is aninteresting research subject. It is important to remember, that climateand other global processes are affected by the conditions of tropicalforests, as the integral units, and not by individual trees or species.
3.5. Additional remarks
The persistent negative intra-seasonal trends coupled withdecreasing SGI and SMI in subsequent years, both being proportionalto the original AIGF, suggest that any reports of apparent inter-annualincreases, decreases, or recoveries of vegetation indexes should beinterpreted with caution. Such observations may be due to differentphenological timing (e.g., early vs. late dry season) of calendar-basedcomparison periods (Souza et al., 2005b; Matricardi et al., 2005).
In this paper, we focused on and established a simple, one-wayconnection between selective logging and phenology. Given recentfindings by Asner et al., 2006, it is reasonable to assume that selectivelogging probabilistically causes future selective logging and deforesta-tion. In this case, the long-term phenological impact of selectivelogging can be exponentially greater. This circumstance magnifies theneed to further evaluate, model, and incorporate the post-loggingphenology of tropical forests into carbon and climate models. In thisway, one will be able to quantify the relationship between phenolo-gical impact and carbon budget at various spatial and temporal scales.
4. Conclusions
This paper presents the first large-scale study on how selectivelogging affects forest phenology, using satellite remote sensing withhigh temporal resolution. Our analyses suggest that selective loggingleading to just 5–10% canopy damage has a significant and long-
lasting effect on forest phenology. Logging progressively reducescanopy moisture content and greenness over the course of the dryseason. Although the overall intra-seasonal impact on greenness andwater content tends to diminish after three years, apparently due tocanopy closure and the reduction of the gap area fraction, there is noevidence that phenology undergoes such a rapid recovery. Observa-tions of changes in forest phenology demonstrate that selectivelogging impacts can be monitored over large spatial extents, apossibility hinted at in field-based studies. Monitoring such changesmay lead to quantifying the impact on carbon flux and sequestrationand determine whether it is of sufficient magnitude to alterbiosphere–atmosphere exchange of CO2, water vapor, and energy inthe logged regions of the Amazon basin. Altered functionality oftropical forests could, in turn, affect a range of biogeochemicalprocesses in the region, and have cascading effects on the regionalclimate system. Our findings highlight the need for future remotesensing and field-based studies to establish the mechanisms andconsequences of these phenomena and to understand their regionaland possibly global ecological importance.
Acknowledgments
Wewish to acknowledgeNASA support under grant #NNG04GK34G,“Studies of biosphere–atmosphere interactionswith a GCMwithMODISSpectral Resolution” and NASA LBA grant NNG06GE32A. We thank theanonymous reviewers whose constructive comments substantiallyimproved the manuscript. Thanks also go to Dr. Michael Keller fromUniversity of New Hampshire for his helpful comments on an earlierversion of the paper. We especially thank Mr. George Scheer and Mr.Lawrence Ross for systems administration and other computationsupport that enabled us to analyze these data.
Appendix I.: Discrete orthogonal polynomials used in theseasonal model
If D+1 points tk divide an interval [a, b] into D=2d equalsubintervals, so that:
tk =a + b
2+ kΔt; where Δt =
b − aD
; k = 0; F1; F2; :::; Fd
and if variable x denotes
x =2t − a + bð Þ
b − ad;
then the orthogonal polynomials en(t) are determined by thefollowing formula:
en tð Þ =Xn
k=0
−1ð Þn + k n + kð Þ 2k½ � d + xð Þ k½ �
k!ð Þ2D k½ � ; n = 0; :::;D
where
z k½ � = z z − 1ð Þ z − 2ð Þ::: z − k + 1ð Þ; z 0½ � = 1 zz0ð Þ; 0 k½ � = 1 for k = 1:2:::
For example, first two orthogonal polynomials are given by
e0 tð Þ = 1; e1 tð Þ = xd;
It is straightforward to show that
XD
k=0
en tkð Þem tkð Þ = 0 for n ≠ m;
which is known as the orthogonality definition.
2440 A. Koltunov et al. / Remote Sensing of Environment 113 (2009) 2431–2440
References
Asner, G. P., Broadbent, E. N., Oliveira, P. J. C., Keller, M., Knapp, D. E., & Silva, J. N. (2006).Condition and fate of logged forests in the Brazilian Amazon. Proceedings of theNational Academy of Sciences of the United States of America, 103(34), 12947−12950.
Asner, G. P., Keller, M., Pereira, R., Zweede, J. C., & Silva, J. N. (2004). Canopy damage andrecovery after selective logging in Amazonia: Field and satellite studies. EcologicalApplications, 14(4), S280−S298.
Asner, G. P., Knapp, D. E., Broadbent, E. N., Oliveira, P. J. C., Keller, M., & Silva, J. N. (2005).Selective logging in the Brazilian Amazon. Science, 310, 480−482.
Cheng, Y. -B., Ustin, S. L., Riano, D., & Vanderbilt, V. C. (2008). Water content estimationfrom hyperspectral images and MODIS indexes in Southeastern Arizona. RemoteSensing of Environment, 112, 363−374.
Cochrane, M. A. (2003). Fire science for rainforests. Nature, 421, 913−919.Feldpausch, T. R., Riha, S. J., Fernandes, E. C. M., & andWandelli, E. V. (2005).
Development of forest structure and leaf area in secondary forests regenerating onabandoned pastures in Central Amazonia. Earth Interactions, 9(6), 1−21.
Gao, B. -C. (1996). NDWI — A normalized difference water index for remote sensing ofvegetation liquid water from space. Remote Sensing of Environment, 58(3), 257−266.
Good, P. I. (2005). Permutation, parametric and bootstrap tests of hypotheses, 3rd editionNew York: Springer.
Herrerias-Diego, Y., Quesada, M., Stoner, K. E., & Lobo, J. A. (2006). Effects of forestfragmentation on phenological patterns and reproductive success of the tropicaldry forest tree Ceiba aesculifolia. Conservation Biology, 20(4), 1111−1120.
Huete, A. R., Didan, K., Miura, T., Rodriguez, E. P., Gao, X., & Ferreira, L. G. (2002).Overview of the radiometric and biophysical performance of the MODIS vegetationindices. Remote Sensing of Environment, 83(1–2), 195−213.
Huete, A. R., Didan, K., Shimabukuro, Y. E., Ratana, P., Saleska, S. R., Hutyra, L. R., Yang,W. Z.,Nemani, R. R., & Myneni, R. (2006). Amazon rainforests green-upwith sunlight in dryseason.Geophysical Research Letters, 33(4) Art. No. L06405.
Huete, A., Justice, C., & Liu, H. (1994). Development of vegetation and soil indices forMODIS-EOS. Remote Sensing of Environment, 49(3), 224−234.
IBGE (2005). IBGE: Sistema IBGE de Recuperacao Automatica: Extracao Vegetal.: InstitutoBrasiliero de Goegrafia e Estatistica (available at: http://www.sidra.ibge.gov.br).
Keller, M., Asner, G. P., Silva, J. N. M., & Palace, M. (2004). Sustainability of selectivelogging of upland forests in the Brazilian Amazon: Carbon budgets and remotesensing as tools for evaluation of logging effects. In D. Zarin, J. Alavalapati, F. E. Putz,& M. Schmink (Eds.),Working forests in the tropics: conservation through sustainablemanagement? New York: Columbia University Press.
Matricardi, E. A. T., Skole, D. L., Cochrane, M. A., Qi, J. G., & Chomentowski, W. (2005).Monitoring selective logging in tropical evergreen forests using Landsat: Multi-temporal regional analyses in Mato Grosso, Brazil. Earth Interactions, 9(25), 1−24.
McNabb, K. L., Miller, M. S., Lockaby, B. G., Stokes, B. J., Clawson, R. G., & Silva, J. N. M.(1997). Selection harvests in Amazonian rainforests: Long-term impacts on soilproperties. Forest Ecology and Management, 93, 153−160.
Monteiro, A. L., Souza, C. M., & Barreto, P. (2004). Detection of logging in Amazoniantransition forests using spectral mixture models. International Journal of RemoteSensing, 24(1), 151−159.
Morton, D. C., DeFries, R. S., Shimabukuro, Y. E., Anderson, L. O., Espirito-Santo, F. D. B.,Hansen, M., & Carroll, M. (2005). Rapid assessment of annual deforestation in theBrazilian Amazon using MODIS data. Earth Interactions, 9(8), 1−21.
Myneni, R. B., Yang,W., Nemani, R. R., Huete, A. R., Dickinson, R. E., Knyazikhin, Y., Didan, K.,Fu, R., Negrn Juarez, R. I., Saatchi, S. S., Hashimoto, H., Ichii, K., Shabanov, N. V., Tan, B.,Ratana, P., Privette, J. L.,Morisette, J. T., Vermote, E. F., Roy, D. P.,Wolfe, R. E., Friedl,M.A.,Running, S. W., Votava, P., El-Saleous, N., Devadiga, S., Su, Y., & Salomonson, V. V.(2007). Large seasonal swings in leaf area of Amazon rainforests. Proceedings of theNational Academy of Sciences of the United States of America, 104(12), 4820−4823.
Nobre, C. A., Sellers, P. J., & Shukla, J. (1991). Amazonian deforestation and regionalclimate change. Journal of Climate, 4(10), 957−988.
Okuda, T., Suzuki, M., Adachi, N., Quah, E. S., Hussein, N. A., & Manokaran, N. (2003).Effect of selective logging on canopy and stand structure and tree speciescomposition in a lowland dipterocarp forest in peninsular Malaysia. Forest Ecologyand Management, 175, 297−320.
Olander, L. O., Bustamante, M. C. C., Asner, G. P., Telles, E., & do Prado, Z. A. (2005).Surface soil changes following selective logging in an Eastern Amazon forest. EarthInteractions, 9(4), 1−19.
Pongratz, J., Bounoua, L., DeFries, R. S., Morton, D. C., & Anderson, L. O. (2006). Theimpact of land cover change on surface energy and water balance in Mato Grosso,Brazil. Earth Interactions, 10(19), 1−17.
Ray, D., Nepstad, D., & Moutinho, P. (2005). Micrometeorological and canopy controls offire susceptibility in a forested Amazon landscape. Ecological Applications, 15(5),1664−1678.
Roujean, J. -L., Leroy, M., & Deschamps, P. Y. (1992). A bi-directional reflectance model ofthe Earth's surface for the correction of remote sensing data. Journal of GeophysicalResearch, 97(D-18), 20455−20468.
Schaaf, C. B., Gao, F., Strahler, A. H., Lucht, W., Li, X., Tsang, T., Strugnell, N. C., Zhang, X.,Jin, Y., Muller, J. -P., Lewis, P., Barnsley, M., Hobson, P., Disney, M., Roberts, G.,Dunderdale, M., Doll, C., d'Entremont, R. P., Hu, B., Liang, S., Privette, J. L., & Roy, D.(2002). First operational BRDF, albedo nadir reflectance products from MODIS.Remote Sensing of Environment, 83(1), 135−148.
Schulze, M., & Zweede, J. (2006). Canopy dynamics in unlogged and logged forest standsin the eastern Amazon. Forest Ecology and Management, 236(1), 56−64.
Souza, C., & Baretto, P. (2000). An alternative approach for detecting and monitoringselectively logged forests in the Amazon. International Journal of Remote Sensing, 21(1), 173−179.
Souza, C., Firestone, L., Silva, L. M., & Roberts, D. (2003). Mapping forest degradation inthe Eastern Amazon from SPOT 4 through spectral mixture models. Remote Sensingof Environment, 87(4), 494−506.
Souza, C., & Roberts, D. (2005). Mapping forest degradation in the Amazon region withIkonos images. International Journal of Remote Sensing, 26(3), 425−429.
Souza, C., Roberts, D., & Cochrane, M. A. (2005). Combining spectral and spatialinformation to map canopy damage from selective logging and forest fires. RemoteSensing of Environment, 98(2–3), 329−343.
Souza, C. M., Roberts, D. A., &Monteiro, A. L. (2005). Multitemporal analysis of degradedforests in the southern Brazilian Amazon. Earth Interactions, 9(19), 1−25.
Stone, T. A., & Lefebvre, P. (2000). Using multi-temporal satellite data to evaluateselective logging in Para, Brazil. International Journal of Remote Sensing, 19(13),2517−2526.
Strahler, A. H., Lucht, W., Schaaf1, C. B., Tsang, T., Gao, F., Li, X., Muller, J. -P., Lewis, P., &Barnsley, M. J. (1999). MODIS BRDF/Albedo Product: Algorithm theoretical basisdocument. version 5.0. available at: http://modis.gsfc.nasa.gov/data/atbd/atbdmod09.pdf.
Ter Steege, H., Boot, R., Brouwer, L., Hammond, D., Van der Hout, P., Jetten, V. G., Khan, Z.,Polak, A. M., Raaimakers, D., & Zagt, R. (1995). Basic and applied research for soundrain forest management in Guyana. Ecological Applications, 5, 904−910.
Tucker, C. J. (1979). Red and photographic infrared linear combinations for monitoringvegetation. Remote Sensing of Environment, 8, 127−150.
White, L. J. T. (1994). The effects of commercial mechanised selective logging on atransect in lowland rainforest in the Lope Reserve, Gabon. Journal of TropicalEcology, 10(3), 313−322.
Xaio, X. M., Hagen, S., Zhang, Q. Y., Keller, M., & Moore, B., III (2006). Detecting leafphenology of seasonally moist tropical forests in South America with multi-temporal MODIS images. Remote Sensing of Environment, 103(4), 465−473.
Xaio, X. M., Zhang, Q. Y., Keller, M., Saleska, S., Hutyra, L., De Camargo, P., Wofsy, S.,Frolking, S., Boles, S., Keller, M., & Moore, B., III (2005). Satellite-based modeling ofgross primary production in a seasonally moist tropical evergreen forest. RemoteSensing of Environment, 94(1), 105−122.
Young, T. P., & Hubbell, S. P. (1991). Crown asymmetry, treefalls, and repeat disturbanceof broad-leaved forest gaps. Ecology, 72(4), 1464−1471.
Glossary
AIGF: area-integrated gap fractionAIGF0: area-integrated gap fraction observed in 2000BRDF: bi-directional reflectance distribution functionETM+: Enhanced Thematic Mapper PlusEVI: Enhanced Vegetation IndexINPE: Instituto Nacional de Pesquisas Espaciais (Brazilian National Institute for
Space Research)MODIS: Moderate Resolution Imaging SpectroradiometerNBAR: Nadir BRDF Adjusted ReflectanceNDVI: Normalized Difference Vegetation IndexNDWI: Normalized Difference Water IndexPRODES: Program for Monitoring Deforestation in the Brazilian AmazonRMSE: root mean squared errorSGI: Effective Intra-Seasonal Greenness ImpactSGIT: Effective Intra-Seasonal Greenness Impact TrendSMI: Effective Intra-Seasonal Moisture ImpactSMIT: Effective Intra-Seasonal Moisture Impact TrendVI: vegetation index
