Remote Sensing Multisensor Ndvi-based Monitoring of the Tundra-taiga Interface (mealy Mountains, Labrador, Canada)

The analysis of a series of five normalized difference vegetation index (NDVI) images produced information about a Labrador (Canada) portion of the tundra-taiga interface. The twenty-five year observation period ranges from 1983 to 2008. The series composed of Landsat, SPOT and ASTER images, provided insight into regional scale characteristics of the tundra-taiga interface that is usually monitored from coarse resolution images. The image set was analyzed by considering an ordinal classification of the NDVI to account for the cumulative effect of differences of near-infrared spectral resolutions, the temperature anomalies, and atmospheric conditions. An increasing trend of the median values in the low, intermediate and high NDVI classes is clearly marked while accounting for variations attributed to cross-sensor radiometry, phenology and atmospheric disturbances. An encroachment of the forest on the tundra for the whole study area was estimated at 0 to 60 m, depending on the period of observation, as calculated by the difference between the median retreat and advance of an estimated location of the tree line. In small sections, advances and retreats of up to 320 m are reported for the most recent four-and seven-year periods of observations.


Background and Rationale
Vegetation cover mapping and monitoring are essential components of knowledge acquisition on the impact of climate changes on Arctic tundra ecosystems [1], part of which is the tundra-taiga interface (TTI), or tree line.Satellite and airborne Normalized Difference Vegetation Index (NDVI) images provide green biomass distribution information in proximity to, and north of, the tree line [2][3][4][5][6].Increased NDVI values for the tundra are mainly related to shrub growth especially with low leaf area indices [7] or fractional area coverage below 40% [8].Vegetation indices and thermal images are also valuable inputs to climate change-related atmosphere constituent and air-temperature models [9][10][11][12][13][14][15].
Multitemporal vegetation index images of northern regions can effectively be built in continental scale mosaics [13,[16][17][18][19] or regional scale time series [20].However, baseline vegetation maps must consider variable scales, or spatial resolutions, in order to benefit from multi-disciplinary data integration that helps better understand ecosystems.Continental scale Advanced Very High Resolution Radiometer (AVHRR) image-based NDVI values may be calibrated with other, higher spatial resolution images, such as those obtained from Landsat, the Moderate Resolution Imaging Spectroradiometer (MODIS) and the Satellite pour l'Observation de la Terre (SPOT) [21][22][23].The sensor's matching passbands allow consistent red and near-infrared spectral information to be extracted [19,[24][25][26], to which the NDVI is based [27].
The relationships of vegetation distribution with environmental factors such as soil moisture [5], permafrost depth [3], seasonal temperature [10,13,14,28,29], atmospheric carbon dioxide [9], geomorphology [6] and land cover [30] are investigated from a wide range of spatial resolutions.Generally, for these studies, the vegetation data come from high temporal resolution (daily to weekly) and coarse spatial resolution (1 km to 1 deg.)AVHRR images, while synchronous environmental data are collected from medium to high spatial resolution images, maps or field data.On the other hand, higher spatial resolution images provide information necessary to understand structural components of the TTI, such as tree islands, functional plant type heterogeneity and forest annual encroachment rate, which is of the order of a few to tens of meters [1,[31][32][33].
Regional scale investigations in Northern Canada [30], Northern Alaska [20,28,[34][35][36], Nunavut [22,37], Manitoba [38] and Siberia [39,40] outline the importance of medium (about 30 m) to high (about 4 m) spatial resolution images for vegetation distribution and related climate change studies; while recent work by [41] shows the benefit of cross-sensor data-based monitoring for socio-ecological studies in the Arctic.However, regional scale vegetation distribution and time series are lacking for the TTI [42].

Objective
This paper presents the changes observed in a TTI from a multisensor NDVI series.The spatiotemporal variations are assessed through the comparison of five multispectral medium resolution images recorded between 1983 and 2008 for an area of the Mealy Mountains, Labrador, Canada.The analyses aim at producing NDVI images that show the green biomass distribution and changes, while considering the impact of atmospheric and cross-sensor inconsistencies.In addition, the multitemporal dataset help describe the spatial characteristics of the areas where the NDVI categories stayed the same over time.Finally, the classified NDVI images are used to derive an estimated tree line position.

Study Area and Data
The study area is located in the Mealy Mountains of Labrador, Canada, approximately 30 km off the south shore of Lake Melville and 100 km east-northeast of Happy Valley-Goose Bay (Figure 1).The study area covered by the image series is square-shaped, centered at 53.60°N and 58.  ) cover most of the understory and forest gaps.Dwarf birch, sparse tree islands and krummoltz patches grow at higher elevations.Lichens (Cladonia spp., Cladina spp., and Cetraria spp.) are found throughout the forest understory, as well as at higher elevations where low lying evergreen shrubs and rocky ground also occur [46,47].
The TTI is interchangeably named tundra-taiga boundary, ecotone, zone or interface [48,49].For this study, the TTI is defined as the transition zone between tundra-dominated and taiga-dominated land covers.
The study area is of interest to the inter-disciplinary project 'Present processes, past changes, spatiotemporal dynamics (PPS) Arctic Canada', an International Polar Year initiative.In this context, the aim of the remote sensing application is to contribute to a research effort in studying the effects of climate change on the tree line at Arctic sites around the globe.
Field data were collected in July, 2008 for the research program on the 'Impact of climate change on alpine habitats and treeline' [50], which is part of the PPS Arctic Canada initiative.The protocol for low Arctic or semi-Arctic tundra sites [1] was applied where data were collected in a cluster of five 1-m 2 quadrats distributed within 63 randomly-designated 30 × 30 m 2 sampling areas.Four of the five quadrats are located in each cardinal direction at a distance of 10 m from the center of the sampling areas, while the fifth plot is randomly selected [1].The percent cover for the five quadrats was compiled to reveal a dominating land cover in each sampled area, where photosynthesis producing tree species (spruce, fir, alder or birch) at least 2.5-m tall were interpreted as taiga, while moss, shrubs or rock were categorized as tundra.The UTM Zone 21-North geographic coordinates of the center of the 30 × 30 m 2 areas allowed locating them on the georeferenced remote sensing images (see Section 6.1.1).
The remote sensing image collection is aimed to include matching medium spatial resolution data from the red and near-infrared spectral bands, from the longest possible multisensor time series, recorded during the vegetation-growing season.These bands are essential to capture the vegetation spectral signature required for computing the NDVI for each year (t) (Equation ( 1)) [51,52].Also, the image selection was based on the combination of best possible synchronicity, or near anniversary date, and lowest cloud cover as catalogued by [53,54] and as inspected visually, respectively.NDVI t = (near infrared t − red t )/(near infrared t + red t ) The series includes five images recorded from different earth observation satellite platforms over a 25-year period.The acquisition of the Landsat (L)-4 Multispectral Scanner (MSS), L-5 Thematic Mapper (TM) and L-7 Enhanced TM+ (ETM+) images of 1983, 1992 and 2001, respectively, was followed by a 2005 Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) image and a 2008 SPOT-5 High Resolution Visible and Infrared (HRVIR) image.The different sensors were used because of a lack of imagery from common sensors.From the earliest image available for the study area, which was recorded in 1983, the image set was gathered to represent each decade, until 2008, when field data was collected.The ETM+ orbital path missed covering of a small portion, approximately 10%, in the northwest corner of the study area where some of the highest terrain elevations occur.
All images but the MSS, which is 6-bit, are in 8-bit quantization level.The spatial resolution is another varying figure of merit across the multisensor image set as it ranges from 80 m with MSS to 10 m with HRVIR (Table 1) [26,55,56].

NDVI Interpretation
The forest canopy fraction and vegetation chlorophyll content can have a significant impact on the NDVI values recorded for the TTI.The forest spectral response is a result of a combination of the reflectance and transmission of the tree canopy and the background reflectance [57].However, the differing vegetation chlorophyll content of black spruce and larch tree species, as well as that of moss and sedge ground cover species cause the NDVI values not to be solely dependent on canopy greenness [57].Spectral signatures recorded for the TTI are for land covers that may range from treeless tundra, including exposed rock, to an open forest cover, the taiga, with a canopy closure of about 25%; with this transition, the proportional composition of herbs, lichen, and bare soil is reduced to the benefit of trees and shrubs [23].
Larch, dwarf birch and graminoid species spread out into the tundra as a response to climate warming at northern latitudes and higher terrain elevations [32,58].The remote sensing derived NDVI is an effective estimation of larch productivity [59] and leaf area index [60].It was shown by [61] that most of the variance in TM and HRV simulated NDVI values for a region of Alaska is the effect of the photosynthetic biomass fraction in tundra tussock.The presence of deciduous shrub foliage such as willow and birch, or live graminoids and forbs considerably increase NDVI values [20,22,35].These can rise to above 0.75 for tundra with shrub, or other photosynthetic biomass content, while they typically range between 0.50 and 0.75 for mossy tussock, non-shrub tundra ensembles' [6,35,40,61].Also, [62] report peak NDVI values on the order of 0.45 to 0.55 for Alaska's moist acidic tundra (subzone E).Evergreen shrub phytomass is not significantly correlated with NDVI values [35], neither are arctic landscapes dominated by moss, lichen and other low-or non-photosynthetic biomass [35,38,61].Conversely, [63] and [64] report visible and near-infrared reflectance signatures for low-chlorophyll moss and lichens that are influenced by the chlorophyll content and other pigmentation.

Expected Cross-Sensor NDVI Differences
Three main issues are expected to affect the ability to monitor NDVI-related vegetation changes over the period of observation.First, the spectral passband differences across sensors, second, the seasonal temperature cycle, and third, the atmospheric conditions.

Across Sensor Spectral Passbands
The red and near-infrared spectral bands are similar across all sensors except for those of the L-4 MSS, which are wider.Particularly, the MSS's near-infrared sensitivity extends into the visible band and thus includes the vegetation 'red edge.'This implies that lower reflectance values may yield lower NDVI values.The TM, ETM + , ASTER and HRVIR infrared bands were designed to only record the vegetation's high reflectance values [17,65].In their work, [24] further examined the relationships between AVHRR and MODIS spectral responses for green vegetation.They found that for high NDVI values, which represent higher green biomass, the ratio calculated from the AVHRR on NOAA-14 (720 to 1,000 nm near-infrared band) could be 7% lower as compared to that from MODIS (841 to 876 nm near-infrared band).NDVI calculated from TM, ETM + , and AVHRR may differ by up to 2%, while these can depart from the MSS-derived NDVI by 5% [66].

Seasonal Temperature-Related NDVI Phenology
The MSS image was acquired the earliest in the growing season and the TM image, the latest.All five images were acquired at least five weeks after the onset of the growing season as indicated by the date when the sum of degree-days above 0 °C reached 300 °C [32] (Table 2).The sum of degree-days is the highest for the ETM+.However, the NDVI phenology is more favorable for the HRVIR image to yield values higher than those from the TM, ETM+, and ASTER images, which were recorded at later dates in the season when the normal temperatures are below 10 °C (Figure 2).The MSS image is unique in that it was acquired during the period of summer temperature maxima.The highest seasonal NDVI values are expected about 10 days to two weeks after the peak seasonal temperature had occurred, and then a gradual decrease of the vegetation index follows, consequent with the temperature change [13,67,68].As a result, the green biomass may not have reached its full potential at the time the MSS image was recorded (see 'MSS' on Figure 2); nevertheless it may have been fuller than when the other images were recorded.The seasonal phenology would have caused the MSS derived NDVI values to be higher than those from HRVIR, followed by those from the other three September images.An evocative approach to assess the NDVI data continuity is by examining temperature anomalies (daily-to-normal temperature difference) on short time periods (of 10 to 40 days) that precede the image acquisition [13].For all five years, temperature anomalies are of the order of ±1 to 2 °C, over a 40-day period before the image acquisition dates (Table 3).A unique situation among the dataset is that over a 10-day period preceding the HRVIR, the temperature anomaly reached +4.1 °C.The smaller temperature anomalies may have caused NDVI variations of about 5%, most likely positive, for the tree line transition location of the study area, while the large short-term temperature anomaly could have caused a seasonal NDVI fluctuation of up to 10% [13].

Atmospheric Conditions
The cloud cover and cloud shadow were identified and masked following the procedure described in Section 6.1.2.The clouds caused 46% and 13% of the study area imaged by TM and ASTER, respectively, to be lost.The MSS and ETM+ are cloud free, while clouds cover only 6% of the HRVIR image (Figure 3).Cumulatively, clouds and their shadows obscure 53% (46% over land, 7% over water) of the study area.The NDVI values from ASTER may have been more affected than the four other images by atmospheric water vapor absorption in the near-infrared band.The atmospheric humidity recorded from the Goose A weather station [69] at the time the images were acquired range between 33% and 43% for all images, except for the ASTER, which was recorded as the atmospheric humidity reached 73%.The effect of water vapor is to reduce the NDVI value, especially in the higher end of its dynamic range [24].
The across-image NDVI differences due to variations of the atmospheric optical thickness are considered negligible as the weather station reports equal values of visibility during all image acquisitions.However, the visibility that was recorded is 24.1 km, approximately corresponds to an aerosol optical thickness of 0.2 of clear sky [70].This would have caused the NDVI values above about 0.4 to be underestimated (see Figure 5 in [70]).

Image Analysis
The image analysis included preprocessing and extraction of information about the TTI.The preprocessing considered first, a geometric correction that brought all five images to a common reference system, spatial resolution and study area coverage.Second, cloud and water masks were produced, which served as constraint layers for only the cloud-free land portion to be considered in further analyses.Third, atmospheric corrections and cross-sensor calibrations were intended to generate a consistent image dataset.
Information extraction consisted of computing and classifying the NDVI images.Then, the spatial structure of a three-class scheme was examined across the multitemporal dataset.Finally, a linear-form representation of the TTI was derived and its change of position was assessed.

Geometric Corrections
An image-to-image geometric correction of the 1983, 1992, 2005 and 2008 data to the 2001-ETM+ orthoimage (UTM Zone 21-North and North American Datum (NAD) 1983) took the multitemporal set to a 30 m pixel size, which complies with the current remote sensing image analysis protocol for mapping Canada's arctic vegetation [1].
The TM, ASTER and HRVIR images' initial spatial resolutions are the same or finer than that of the reference ETM + image (Table 1).The nearest neighbour interpolation was applied for those images with which the ground control points could be located based on visual interpretation of similar brightness values in corresponding spectral bands.However, the MSS image spatial resolution was much coarser than that of the output matrix.In order to produce a gradual transition from the input to the output resolution, the path-projected MSS pixels were expanded upon by a factor of four in row and column directions, to a 20 m pixel size, and the bilinear interpolation resampling function was applied.The pixel replication facilitated the relative location (i.e., 'middle' or 'edges') of ground control points within the same-brightness values areas corresponding to the initial 80 m image units.
The number of ground control points for each image is 16 for the 1992, 2005 and 2008 images and 12 for the 1983 image.GCPs included small islands and lake shorelines with unique profiles, as anthropogenic features such as road intersections, were not present in any of the images.The GCPs were uniformly selected throughout each image.Each resulting geo-corrected image has a total root-mean-square error equivalent to less than a half of their initial spatial resolution (37, 15, 3, and 7 m for MSS, TM, ASTER and HRVIR, respectively).

Cloud, Cloud Shadow, and Water Masks
The 1983 MSS image served as the cloud-free image standard for the discrimination of clouds and cloud shadows on the 1992, 2005 and 2008 images.The image ratio change detection analysis [71,72] consisted of applying a two-date red spectral band image ratio of the cloud-free image to the cloud-affected image.The overlay operation enhanced the clouds and their shadow.On the cloud-enhanced image histogram, the lowest and highest values characterize the cloud shadows and clouds, respectively.Once classified, these areas were expanded upon with a five-pixel buffer.This buffer width was chosen because it included thinner clouds that were identified by visual inspection of the images, but that had not been included in the automatically classified clouds category.
NDVI images are relevant only to the land portion of the image, therefore justifying the application of a water mask.The near-infrared band bimodal histogram, which displays distinctive distributions for water and land [73], is classified with reference to the threshold reflectance value that recorded the lowest histogram frequency between the water and land values.
The clouds, cloud shadows and water mask were integrated through image overlay to each image of the set and, separately, a cumulative 'no data' mask was applied to the images.

Atmospheric Corrections and Radiometric Cross-Calibrations
Continuity of NDVI values ideally relies on sensor pairwise cross-calibration of the individual bands reflectance values, which may vary in relation to the land cover [25,65].The radiometric preprocessing included, first, a calibration of the individual images and spectral bands to the 'at-satellite apparent reflectance' value, with the integration of a Rayleigh scattering path radiance equivalent (i.e., dark-pixel subtraction) correction.Second, a multitemporal cross-calibration adjusted the reflectance values of corresponding spectral bands across the five images.Both calibrations methods were used in order to minimize the potential errors due, on the one hand, to atmospheric moisture conditions that affected the images and that would not have been all accounted for in a path radiance equivalent correction, and on the other hand, the difficulty to find a calibration points (same land cover) in images that were collected over a period of 25 years and varying cloud covers.
The radiance calibration values for the three Landsat images were documented by [26] while the HRVIR parameters were contained in the image metadata.No calibration parameters accompanied the ASTER image, therefore, in this case, the individual spectral band calibration was omitted and the raw image was input to the cross-calibration.The individual spectral band calibration incorporated atmospheric corrections based on the dark object subtraction method to minimize the effect of Rayleigh scattering [74].The image histogram non-zero minimum brightness values (BV) provided the haze equivalent estimates.
The multitemporal relative correction, or cross-calibration, aimed to normalize to a reference image, the at-satellite apparent reflectance of its corresponding spectral band.The 2001-ETM + image served as the reference image due to the cloud-free conditions during which it was acquired.In addition, the ETM + red and near-infrared passbands are best matched to the corresponding bands of the other images in the dataset.The calibration is based on regression equations obtained from 31 points, each chosen to represent a pseudo-invariant feature (PIF) [75] throughout the study area.PIFs for this particular area included barren surfaces and deep water bodies, as these features are shown to have little variation in spectral reflectance across the multitemporal dataset.The spectral characteristics of a PIF must be consistent between the reference image and the image being calibrated.An example of the linear regression graph between the ETM+ image and the SPOT image is presented in Figure 4.The NDVI computed from the near-infrared and red calibrated images of each year (Equation ( 1)), provided a representation of the green vegetation cover distribution.A red reflectance decrease and an infrared reflectance increase translate into an increased NDVI.
The NDVI images were classified into ordinal categories that represented the vegetation cover while incorporating the anomalies identified to the sensors' spectrometry, atmospheric conditions and to the vegetation phenology throughout the multitemporal dataset.The areas for which data are cumulatively available on all five images (cloud-, cloud shadow-free and in-image path) were used to set the NDVI classification boundaries.The positive NDVI values were grouped into three categories based on quartile distributions.The interquartile values (middle 50%) of each year-distribution were grouped in the class 'intermediate'.The values included in the lowest and highest quartile ranges were assigned to the 'low' and 'high' NDVI classes, respectively.
The ecological significance of the NDVI classes was evaluated using the Fisher's exact probability test [76].The null hypothesis is of no difference in the 'low' and 'high' NDVI categories considering the field data representing 'tundra' and 'taiga' land cover types as defined in Section 3 for which field data were available.

NDVI Spatial Distribution and Within Class Variation
Information about the horizontal distribution of the TTI was extracted by overlay of the five image set to identify the areas that are, on at least four of the five images, repeatedly in one of the three NDVI classes.This analysis caused patterns of consistent low, intermediate and high vegetation indices to emerge.
The within-class median NDVI value for each of the three consistent NDVI categories was computed.The results of this analysis are represented and put in context of the errors that were attributed to spectral resolution and short-term temperature anomalies.

Estimating the Location of the Tree Line
Tentatively, the intermediate NDVI class represents a transition class between the taiga and the tundra land cover classes.This class is expressed as the TTI, i.e., where the tundra and taiga environment interact.The transition from taiga to tundra is gradual because it occurs along a sizable distance that forms a corridor-shaped area.The TTI, as represented by the intermediate NDVI class, cannot be monitored based on the information that we have extracted and interpreted from the dataset because it is spatially fragmented (Table 5).For example, the maximum area size for continuous intermediate NDVI patches is 7 ha, with a median size of 1 ha which correspond to only a few TM or MSS pixels.In addition, the field data were collection protocol focused land cover description within 30 × 30 m 2 sampling areas.Therefore, it is practical to further this study by estimating a 'tree line' halfway between the low NDVI and high NDVI classes which were validated using field data collected for tundra and taiga land cover types, respectively.This tree line was estimated to help assess the magnitude of changes.The tree line estimation process flowchart is presented on Figure 5.The area perimeter of low and high NDVI served as anchor to calculate the radiating Euclidean distance and to automatically identify the pixels that are halfway between them, at which location an 'estimated' tree line is drawn.The halfway distance was identified by subtraction overlay where the output image value of zero was assigned to pixels situated where the distance calculated from the low NDVI is equal to the distance from the high NDVI.With each image set, the distance was automatically calculated from the tree line on the earliest image date (1983 or 2001) to the displaced tree line segments on the second and third images (2001 and 2008, or 2005 and 2008).The displacement distances were calculated by overlay of the estimated tree line-migrated segment image and a distance gradient from the earliest date's estimated tree line image.

Results
Table 6 presents the descriptive statistics of the five NDVI images.Relatively low mean and maximum NDVI values emerged from the 1983 MSS and the 1992 TM data.The NDVI values calculated from the 2008 HRVIR image are the highest, with an average value of 0.60.It is about 0.05 to 0.20 higher than that of the average NDVI computed from the other images.
The quartile distribution-based NDVI classes are listed in Table 7 and the most consistently classified areas are displayed on Figure 6; a consistent area is where same one of the three NDVI class is present on at least four of the five images.Areas where the low and high NDVI remained similar over the 25-year period have the largest and most variable unit sizes (Table 5).Conversely, the intermediate NDVI category is highly fragmented into parcels that are of the order of 1 to 3 ha in size.The within-class median NDVI value computed for each of the three consistent NDVI categories show an increasing trend from 1983 to 2008 (Figure 7).During this period, the median values for the low, intermediate and high NDVI classes increased by 42%, 39%, and 34%, respectively.These changes are robust to the estimated errors due to a wider passband on the 1983 image (7%) and short-term temperature anomalies (5% for all images and 10% for the 2008 image).
Table 8 presents the contingency tables and the Fisher's exact test results.The null hypothesis was rejected for all five images, meaning a significant association (α = 0.05) between the classes that describe the lowest and highest NDVI values and the field-observed vegetation cover.Given the NDVI values involved that were documented in the literature (see Section 4) [6,35,40,61,62] and the Fisher's test results, the low, intermediate, and high NDVI classes were matched with the tundra, TTI, and taiga land cover types, respectively.The multitemporal image on Figure 8 illustrates the 1983 to 2008 changes.Lines of different colours represent the tree line location at different dates and the black line indicates that no change occurred over the period of observation.For visual reference, the green-gridded feature in the southeast portion of the image is the area of high NDVI for 1983, which is the initial date of the three-image series.
Segments of a few kilometers in length of the estimated tree line keep a steady position with changes that do not exceed 90 m over the period of observation.At other locations it is unstable as it alternatively records advances and retreats, as well as displaying steady advances or retreats.At locations in the southeast portion, in the low lands of the study area, the estimated tree line makes advances of the order of 200 to 300 m when comparing the 1983 and the 2008 images.Also in this area, the estimated 1983-tree line showed several discontinuities (evidenced by a multitude of circling tree line segments) that were filled-in as represented on the 2001 and 2008 images.
The multitemporal image on  Changes of the tree line estimated at halfway between the low and high NDVI classes, indicate a net progression of 10 to 60 m over the different periods of time.Locally, the variations are much larger than the image's spatial resolution.The tree line segments that are subject to gradual migration or retreat depart by up to 320 m from the estimated position of the 1983 or 2001 tree line.The most stable portions of the tree line, which have not shown change of more than 90 m, are of the order of 3 to 4 km in length (Figures 8 and 9).
The clouds are not shown on the composite images because the tree line location is estimated from the low and high NDVI class patches, through this process any other class would not be considered.The 1983 image is cloud free, as well as the 1992 image.The clouds on the 2008 image do not interfere because they are in the southeast corner of the study area (Figure 3) and away from the transition zone.Small clouds may have affected the 2005 estimated tree line on the short-term change image by spatially interrupting the distribution and producing smaller patches of low and high NDVI from which the tree line was calculated.

Discussion
Changes in green biomass, represented through the NDVI, and its spatial distribution were mapped from a multitemporal image series acquired by five different sensors.The analysis resulted with ordinal NDVI categories that incorporate the variations due to acquisition asynchrony caused by the seasonal temperature anomalies and cross-sensor spectral inconsistencies.The quartile-based NDVI classes are representative of land covers found near the tundra-taiga interface.Throughout the 25-year period of observation, all three classes display an increase of the trend in NDVI values.These are of the order of 40%.In addition, an estimated tree line was derived, which allowed to illustrate and quantify long-and short-term changes in the landscape.Some of the sources of NDVI anomalies reviewed in Section 5 may have affected the dataset.In particular, lower NDVI values in the 1983 image may be attributed to the different infrared passband on the MSS sensor.The TM-derived NDVI image displays the smallest coefficient of variation as well as the lowest average value.These may have been caused by a lower level of incident radiation due to extensive cloud coverage during the 1992 image acquisition.Yet, the low and high NDVI categories were tied with the field-recorded land cover descriptions.
The highest NDVI values calculated from the 2008 HRVIR may have been caused by the short-term temperature anomalies prior to the image acquisition date.Nevertheless, the phenological state at the time of the MSS image was recorded, one month earlier in mid-July, also works in favor of a greener land cover.This might explain on the long-term change image that some portions of the estimated tree line have receded.
An important issue of multitemporal image analysis is that the validation of historical images, which are these images that were recorded before one engages in a monitoring study, have no coincident field data or very little if fortuitously it was collected for other purposes.The Fisher's test results reported in Table 8 are all significant and there is no apparent unbalance of the distribution of the observations in the contingency tables for the different images, including that of coincident HRVIR image.This is an indication that the field data were collected in area where no important change had occurred.In addition, the contingency tables reveal a mix of tundra and taiga field observations in the intermediate NDVI category; because of this it is anticipated that this class gives a correct representation of the TTI on the image.Noticeably, this area contains a larger proportion of tundra samples in 1983 (14 of 22 observations, or 64%) and 1992 (61%), while a more even representation of tundra and taiga field records are found in the intermediate NDVI class on the 2001, 2005 and 2008 images.
The categories of low and high NDVI occur in patches that have dimensions that are much larger than the dataset's coarsest resolution (that of MSS with 80 × 80 m 2 or 0.64 ha pixel size).Also, more than half the parcels forming the transition zone (intermediate NDVI class) have sizes of about 1 ha.Therefore, as expected, the MSS spatial resolution is inappropriate for monitoring the progression of features such as tree islands and encroachments that are found at the boundary and within the TTI.For this purpose, neither is appropriate the TM image 30 m resolution, with about 11 pixels per ha.
The image analysis revealed increase of the TTI density.The 1983-to-2008 change image (Figure 8) indicates in-filling of areas that were not in the high NDVI class in 1983.The in-filled areas are shown by the 1983-line (green) forming closed loops that are not duplicated by the other tree lines (yellow and red).These patterns take place on the long-term change image only.A similar dynamics, whereby changes occurred from sparse to open, and from open to normal stand was reported by [32] for the period of 1980 to 2000 in forest-tundra larch forests Ary-Mas, Russia.
The investigation of the TTI as an area feature that corresponds to the intermediate NDVI category must be based on spatial and temporal resolution images that have not been available until just recently.Given the 1 ha area formed by intermediate NDVI patches, and considering a 3 × 3 pixel array for being able to observe a 1-pixel migration, 30 m is the coarsest spatial resolution appropriate for detecting variations within the transition area.In light of this evaluation, the spatial resolution prescribed in the protocol for mapping Canada's arctic vegetation [1] may have to be reconsidered when focusing on regional scale observations of the tundra-taiga interface.
An uneven advance of the tree line is observed from the result of this study.While several segments may not have changed at all over the years, some portions have retreated and others display 270 m (68 m per year from 2001 to 2005) to 320 m (46 m per year from 2001 to 2008) horizontal progression in the direction of the tundra.These values are lower than the rates of 0.2 to 0.4 km per year stated by [48] for Eastern Canada.
The difference between the advance and retreat median values concludes to a net advance of 9 m per year, for the period from 2001 to 2008 (Table 9: 150 m advance less 90 m retreat over a 7 year period).This is consistent with the advances of 3 to 11 m per year for the period of 1980-2000 in northern Russia [32], but it is much lower than the mean advance of 378 m per year (2.5 km from 2000 to 2009) reported for the pine forest line in northern Norway [77].The variations observed for the Mealy Mountains study are attributed to site conditions, as seedling in the TTI in Labrador can be affected by the soil wetness, wind exposure, temperature [78], seedbed composition [79] and the topography [80].The impact of these factors, individually or in combination, yet has to be investigated in relation to the information that was extracted from the satellite images.
The revisit frequency of visible and near-infrared image systems is in effect reduced due to cloud cover and atmospheric moisture degraded signal.Such issues prevent the acquisition of high temporal resolution image series that would allow for characterization of the vegetation's annual phenology and to differentiate this from long-term changes.The prospect of using radar image-based vegetation mapping [39,81] circumvents this issue by potentially increasing the temporal resolution, while spatial resolutions available are finer than the 30 m used in this study.Further work must investigate the information that can be extracted from radar images, which is not dependent on the vegetation's photosynthetic activity, as the NDVI is, but on the spatial structure of the components of the tundra and taiga environments and of the transition zone separating them.

Conclusion
This paper presents evidence that the green biomass, represented by the Normalized Difference Vegetation Index (NDVI), has increased between 1983 and 2008, and more importantly between 2001 and 2008, in a portion of the Canadian tundra-taiga interface (TTI).The results of this study are based on the analysis of a multisensor NDVI image series that had been recorded using different spectral and spatial resolutions, and under variable atmospheric conditions.A cross-sensor calibration and the classification of the vegetation index in ordinal classes produced a consistent dataset from which a NDVI-inferred biomass overall increasing trend of about 40% was detected.The change in biomass is more accentuated for the low NDVI than for the high NDVI category, which rendered an increasing trend of 42% and 34%, respectively.
Spatially, a migration of the taiga environment, represented by high NDVI values, toward the tundra occurred.This observation emerged from changes in the location of a tree line derived from the classified images.Overall, no significant difference was detected between the tree line estimated from the 1983 and the 2008 images, while a net advance of 60 m, or 9 m per year, was observed between 2001 and 2008.Horizontal migrations of 130 to 320 m were recorded for short segments of the estimated tree line during the various periods of observations.
In characterizing of the TTI as an area, rather than as a linear feature, the research highlighted that it has a large spatial variability, represented by the patch size, of the low and high NDVI classes.Consequently, the monitoring of the TTI would preferably rely on image spatial resolutions of 20 m (25 pixels per ha) or finer.
A multisensor image set can successfully be used to monitor the TTI and potentially expand to include the latest medium to high spatial resolution multispectral images such as acquired from the Landsat-8, QuickBird, and WorldView-2 platforms.As well, further work for the study area of the Mealy Mountains, Labrador, addresses the transition from multispectral to radar image for monitoring the tundra-taiga interface.
83°W and is bound by the (373,440 mE, 5,935,343 mN) to (383,670 mE, 5,945,573 mN) coordinates of the Universal Transverse Mercator (UTM) Zone 21-North.

Figure 1 .
Figure 1.Location map of the study site in the Mealy Mountains, Labrador, Canada.

Figure 2 .
Figure 2. Image acquisition dates plotted on the climatic average daily temperatures recorded at the Goose A weather station.

Figure 4 .
Figure 4. Linear graph of the red and infrared reflectance of the pseudo-invariant features.This example is for the HRVIR to ETM + calibration.

Figure 5 .
Figure 5. Process flowchart used to extract the estimated tree line as a linear form feature.

Figure 6 .
Figure 6.Areas classified as low, intermediate and high NDVI on at least four of the five images years.The class 'variable NDVI' is for areas that were not consistently (on less than four images) assigned to a same class from image to image.

Figure 7 .
Figure 7. Median NDVI value in the five quartile-based classifications.The values reported are for the area of matching classes on at least four of the five images of the time series.The error bars express the estimated variations due to the spectral resolution and short-term temperature anomalies.A grey line marks the linear trend for each class.

Figure 8 .
Figure 8.Long-term multitemporal composite of the tree line estimated from the 1983, 2001 and 2008 images, with the high 1983 NDVI overlaid as green grid.
Figure 9 illustrates the 2001 to 2008 changes.It is constructed and interpreted in the same manner as Figure 8 except for that the green gridded area represents the 2001 high NDVI class.As observed on the long-term multitemporal image, segments of few kilometers of the estimated tree line have remained stable between 2001 and 2008.In the center area of the image, noticeable advances from high to low NDVI are recorded, especially for the last portion of the time series, between 2005 and 2008.

Figure 9 .Table 9 .
Figure 9. Short-term multitemporal composite of the tree line estimated from the 2001, 2005 and 2008 images, with the high 2001 NDVI overlaid as green grid.

Table 1 .
Red and near-infrared spectral band ranges, and nominal spatial resolution of images.

Table 2 .
Date of image acquisition, date when the sum of degree-days above 0 °C (SD 0 °C) reached 300 °C, and SD 0 °C before image acquisition date (BI).

Table 4
reports the cross-calibration equation parameters and coefficient of determination values.

Table 5 .
Area and variability of area of x-and-y connected patches consistently classified as low, intermediate and high NDVI over the 25-year monitoring period.

Table 6 .
Minimum, maximum, average, and standard deviation (SD) of the positive NDVI values for each image.

Table 7 .
Quartile distribution-based low, intermediate, and high NDVI class boundaries for each image.

Table 8 .
Co-occurrence of sampled field observations with image's low, intermediate (Int.) and high NDVI categories, and no data due to cloud and cloud shadow.Fisher's exact probabilities results provided at the bottom right are for 2 × 2 tables of 'tundra' and 'taiga' field data vs. 'low' and 'high' NDVI image classes, respectively.