OPEN ACCESS

Abstract: In this paper we present a framework for modelling potential species distribution (PSD) of balsam fir [bF; Abies balsamea (L.) Mill.] as a function of landscape-level descriptions of: (i) growing degree days (GDD: a temperature related index), (ii) land-surface wetness, (iii) incident photosynthetically active radiation (PAR), and (iv) tree habitat suitability. GDD and land-surface wetness are derived primarily from remote sensing data acquired with the Moderate Resolution Imaging Spectroradiometer (MODIS) sensor on the Terra satellite. PAR is calculated with an existing spatial model of solar radiation. Raster-based calculations of habitat suitability and PSD are obtained by multiplying normalized values of species environmental-response functions (one for each environmental variable) parameterized for balsam fir. As a demonstration of the procedure, we apply the calculations to a high bF-content area in northwest New Brunswick, Canada, at 250-m resolution. Location of medium-to-high habitat suitability values (i.e., >0.50) and actual forests, with >50% bF, matched on average 92% of the time.


Introduction
Understanding tree species distribution is essential for the sustainable management of forest landscapes and conservation of biodiversity [1][2][3].Tree species distribution, occurrence, abundance, and growth are all affected by the same primary factors of temperature, soil water availability, incident solar radiation, and soil geochemical properties [4][5][6][7][8].Sudden changes in any one of these factors can produce important and unalterable changes in the composition and ecology of landscapes.Specific variables with potential to control actual species distribution include: (i) accumulated growing degree days, (ii) precipitation and redistribution of rainwater or snowmelt as surface or sub-surface flow, (iii) surface desertification (especially in arid to semi-arid regions of the world), (iv) minimum winter temperatures, (v) winter thaw-freeze severity [9], (vi) number of frost-free days, (vii) snow accumulation, (viii) soil acidification, (ix) forest floor thickness, (x) soil composition and compaction, (xi) intra-and inter-species competition, (xii) disturbance regimes [10], and (xiii) time since disturbance.
In order to model tree species distribution, a common approach is to relate biologically-relevant environmental variables to species habitat preferences in GIS [3,5,[11][12][13].However, estimation of biophysical properties for large areas has often been based on the collection of point data and numerical interpolation.In general, the quality of interpolation reflects the quality of the spatial pointdata and the interpolation method used.Many interpolation methods (e.g., kriging, inverse distance weighing, natural neighbours, nearest neighbours, triangulation, cubic splines [7,14,15]) are incapable of addressing variation in the target measurements as a function of extrinsic factors, such as, for example, time, changes in elevation, slope, direction of the prevailing wind, distance from the coastline, and overall synoptic conditions, other than the internal relationships contained in the data due to autocorrelation.Some interpolation methods like co-kriging, can establish relationships with some external factors, but this is usually limited to a few co-variates.In general, interpolation by most of these methods produces unrealistic values in areas with poor data coverage (with >20-km separation between observation points; [16]) and with strong underlying physical gradients, such as in mountainous terrain [17].One solution to improving spatial representation of biophysical variables of landscapes is by using remote sensing (RS) data.RS data provide a near continuous sampling of the earth's surface at practical spatial resolutions and revisitation intervals.
The use of RS data in species distribution studies is not new.For example, (i) Saatchi et al. [18] employed Moderate Resolution Imaging Spectroradiometer (MODIS)-based estimates of leaf area index and percent tree cover, Quick Scatterometer (QSCAT) backscattering, Shuttle Radar Topography Mission (SRTM)-derived digital elevation model (DEM), and Tropical Rainfall Mapping Mission (TRMM)-based rainfall records to model current tree species distribution in the Amazonian basin; and (ii) Zimmermann et al. [19] combined Landsat ETM+ derived indices (i.e., normalized difference vegetation index (NDVI), green vegetation index, wetness index, and soil brightness index) to model species distribution of rare, early successional, and broadleaf-tree species in Utah, USA.In 2003, Turner et al. [20] predicted the widespread development of mapping techniques of landscapebiophysical conditions with MODIS-sensor acquired data (250-m to 1-km resolution, with a 1 to 2-day revisitation cycle).In previous work we have developed several methods to map the biophysical properties of landscapes relevant to understanding tree species distribution in Atlantic Canada, in particular growing degree days (GDD, a temperature-related index of plant growth; [21]) and temperature-vegetation, land-surface wetness index (TVWI; an indirect measure of soil water content, SWC [22]).
In this paper, we provide a framework for modelling potential species distribution (PSD) of balsam fir [bF; Abies balsamea (L.) Mill] as a function of spatial indices of MODIS-based temperature (by means of GDD) and land-surface moisture (TVWI; [21,22]), and incident photosynthetically active radiation (PAR) simulated with the Landscape Distribution of Soil moisture, Energy, and Temperature (LanDSET) model at 250-m resolution [5,17].No attempt is made to address forest and landscapeforming dynamics associated with inter-and intra-species competition, forest succession, species invasion and extirpation, and disturbance in our definition of PSD.Description of the methods and their application in assessing habitat suitability and PSD are given in Sections 2.2 and 2.3.Spatial indices of habitat suitability and PSD are provided here as values ranging between 0 and 1, where 0 coincides with unfavourable site conditions and, thus, low probability of species occurrence, and 1, superior site conditions.As a demonstration of the method, we model habitat suitability and PSD for a bF-dominated region of northwest New Brunswick, Canada, and compare modelled results with the occurrence of high bF-content stands (with >50% bF, by volume) from an existing digital forest-cover map of the region.In this comparison, high bF-content stands are selected to eliminate multi-species effects and the potential masking of relationships between bF occurrence and biophysical characteristics of the landscape.A comparison between PSD and existing bF cover is meaningful in this study because bF, as an aggressive tree species, regenerates very well on its own on appropriate sites.As a result, bF is normally not planted in New Brunswick, except for the production of Christmas trees and for research purposes (D.E.Swift, personal communication, 2009).

General Description of Study Area
The study area is in the northwest of the Province of New Brunswick (NB), which falls in the Atlantic Maritime ecozone of eastern Canada (Figure 1a).The area is characterized by temperate evergreen-deciduous mix forest-cover types.Geographically, the area is located between 46 o 12'N to 48 o 04'N latitude and 65 o 38'W to 69 o 03'W longitude, and is representative of Highland, Northern Upland, and Central Upland ecoregions of NB (delineated as zones 1a, 1b, 2, and 3a in Figure 1b).Balsam fir is the dominant species, with >50% composition, as identified from aerial photographs [23].The area is characterized by its variable topography, with elevations ranging from 100 m to 840 m above mean sea level (amsl).The cool-moist climate of the area is largely influenced by the area's proximity to the Bay of Chaleur in the northeast, movement of continental air masses from the west, and interior highlands.Mean annual temperature and total precipitation in the area varies from 3.5 o C to 6.5 o C and from 900 mm to 1,500 mm [24]. of high bF-content stands (with >50% bF, by volume) from an existing digital forest-cover map (28120 stands, in all).Forest cover in ecoregion 3a is not available for this study.

Data Requirements and Processing
Data employed in this work included (i) RS-based products of GDD and TVWI and estimates of PAR generated with the LanDSET model for the calculation of habitat suitability (HS) and PSD, and (ii) an existing digital forest-cover map (i.e., Forest Inventory Database of the NB Department of Natural Resources, NBDNR; Figure 1b) for species-response equation parameterization (Section 2.3) and validation of PSD.

Growing Degree Days (GDD)
A GDD map at 250-m resolution was previously developed by employing (i) MODIS-based 8-day composites of surface temperature (T S ; at 1-km resolution) and 16-day composites of enhanced vegetation index (EVI; at 250-m resolution) for the April-October period of 2003-2005, (ii) towerbased half-hourly emitted infrared radiation to estimate canopy T S , and (iii) 30 years (1971-2000 period) average point-estimates of GDD from 101 climate stations from across the Atlantic Maritime ecozone (refer to [21], for greater detail).As an initial step, an empirical relationship was established between mean tower-based estimates of T S for the MODIS-image acquisition period (i.e., time between 10:30 am and 12:00 pm) and the daily mean T S .This relationship was subsequently used to estimate daily mean T S from MODIS-derived values of T S at 1-km resolution for each 8-day period and GDD using the standard definition, i.e.: ( ) where a T is the daily mean air temperature and equal to daily mean T S [25][26][27], T base is the threshold temperature, below which plant growth ceases (approximately, 5 o C for bF; [21]), and i = 1 represents the start-day and i = n, the end-day of the growing season.Figure 2 demonstrates the relationship between a T and daily mean T S over a bF forest near the southern limit of the study area (Nashwaak Lake, NB).Statistical comparison of values yielded an r 2 -value = 0.99 and regression equation slope of 1.01 and intercept of 1.98 K.Note that in order to enhance the spatial resolution of the initial GDD map originally at 1-km resolution to 250-m resolution, a data fusion method was developed relating enhanced GDD to seasonally-averaged EVI [21].Correlation between GDD and EVI was demonstrated to be fairly strong (r 2 = 0.87).As a result, the application of seasonally-averaged EVI's to enhance GDD-map resolution was justified [21].To generate the 30-year average of GDD at 250-m resolution, a constant offset of -511 was applied to reflect the 1971-2000, 30-year normal period [21].collected over a predominantly bF forest near the southern limit of the study area.

Temperature-Vegetation Wetness Index (TVWI)
A TVWI map at 250-m resolution was developed to address the availability of surface (soil) water [22].representation of TVWI across variable terrain required elimination of terrain-induced variations in surface air temperature.To do this, MODIS-based estimates of T S were converted to potential surface temperature (θ S ) at mean sea level by applying basic meteorological principles of adiabatic compression and a 250-m resolution DEM of the study area in the calculation of vertical atmospheric pressure and air temperature variations in an assumed neutrally-stratified atmosphere [22].Estimates of TVWI were obtained by interpreting θ S -NDVI scatterplots generated for each 8-day composite (Figure 3).A long-term average TVWI-map was then produced by averaging all of the individual 8day TVWI images.TVWI values coinciding with 0 represented the dry sites and 1, the wetter sites.
The DEM was generated from 3-arc second SRTM height-data acquired from the National Aeronautics and Space Administration (USA) archives.

Photosynthetically Active Radiation (PAR)
In calculating incoming solar radiation and PAR, terrain attributes of: (i) slope; (ii) aspect; (iii) horizon angle; (iv) view factor; and (v) terrain configuration factor were obtained by processing height data from the DEM [5,17].Incoming solar radiation at the top of the atmosphere was partitioned into its direct and diffused components and each was affected differently as it passed through the assumed cloud-free atmosphere and interacted with the underlying terrain before being summed at the surface [17].Hourly values generated were then added to provide a growing season total.Available PAR was assumed to be 45% of calculated incoming solar radiation [24].Finally, PAR was normalized (nPAR = PAR/PAR max ; where PAR max is the available energy on south-facing slopes) in order that the output from the species-response function for PAR (addressed in Section 2.3) would be constrained between 0 and 1, as would the others.

Modelling of PSD
A conceptual framework for modelling PSD is shown in Figure 4. Development of species-specific environmental response functions for GDD, TVWI, and nPAR followed procedures described in the literature (e.g., [5,28]).TVWI and GDD were both modelled as parabolic functions.A parabolic function for TVWI was legitimate because bF trees grow best in moderately wet soils, and less well in excessively wet or excessively dry soils.We used a linear response function for nPAR.As Barden noted [29], a linear relationship was observed to exist between crop yield and a growing-season integration of PAR in shade-tolerant grass stands.Similar trends were observed in shade-tolerant forests [28].
Equation parameters GDD min , GDD max , TVWI min , and TVWI max are the minimum and maximum values of GDD and TVWI, respectively.For bF, equation parameters were set according to the values, GDD min = 700, GDD max = 2000 (after [30]), TVWI min = 0.12, and TVWI max = 0.52.Parameterization of TVWI min and TVWI max was based on the occurrence of high bF-content stands described in the digital forest-cover database and associated values of TVWI (Figure 1; [22]).
In a manner similar to the environmental limiting-factor approach of Acevedo et al. [31], the compound effect of the biophysical variables on habitat suitability (HS) and PSD (raster-image of HS) was represented as a product of species environmental-response functions, i.e.,

Biophysical Variables
Figure 5 illustrates both the spatial and frequency distributions for GDD, TVWI, and nPAR.GDD values were observed to be between 950-1800 from the coolest to the warmest locations in the landscape (Figure 5a,b).These values are consistent with values reported by Environment Canada for GDD calculated with a base temperature of 5 o C (Canadian Climate Normals for NB, 1971-2000; http://climate.weatheroffice.ec.gc.ca/climate_normals/index_e.html,website last accessed on July 2009).In the case of TVWI (Figure 5c,d), a comparison with LanDSET-calculated SWC [5,17] revealed very similar distributions (r 2 = 0.96; [22]).A comparison between field-based measurements of PAR and modelled PAR (Figures 5e and 5f) at two forest-sites in NB using LanDSET, also provided reasonable agreement (r 2 = 0.96; [24]).

Spatial Distribution of PSD in Comparison to Existing Forest-cover
Figure 6 shows the spatial distribution of modelled HS and PSD for bF in northwestern NB.We summarize PSD patterns in the following way: 1. Polygon A shows an area with low-to-medium HS values (i.e., in between 0.25-0.63; Figure 6a,b).This area is prone to lower PAR, increased TVWI due to physiographic position, and relatively high GDD as a result of the area's proximity to the Bay of Chaleur to the north.This area most likely provides a better growing environment for temperate, hardwood species.This is consistent with species distribution maps produced by Godin and Roberts [23].The few bF stands in the area (Figure 6b) support this interpretation.2. Polygon B shows an area that is partially affected by a rain shadow produced by the Central Uplands of NB (producing lower TVWI); dominant HS values in the area are in the medium-to-high value range (i.e., 0.50-0.75; Figure 6a,b).However, the lack of bF stands in the area (Figure 6b) contradicts what we would expect with high HS.This may be related to the fact that bF in most of these stands were affected by spruce budworm attack during 1970-1990 [32]  An overlay analysis carried out to examine the relationship between high bF-content stands (with >50% bF; Figure 1b) and HS values (Figure 6a) revealed that on average 92% of the stands fell in zones of medium-to-very high HS (>0.50;Table 1).Table 1.A summary of an overlay analysis of high bF-content stands with respect to modelled HS; "n" is the number of stands falling in each bF-content category ("%").

Future Considerations
In order to enhance the current capacities of PSD modelling in northwestern NB, we need to consider the following issues: i. Ecoregion 3a (Figure 1b) has limited forest cover indicated compared to the other regions as most of the land in that region is owned by private companies, limiting the evaluation of HS and PSD for that area.ii.An important site variable, the soil (e.g., [5,12]), is currently not incorporated in the definition of PSD.Existing soil maps are based on geological formations and, as a result, they tend to be unrepresentative at smaller spatial resolutions.iii.Implementation of the framework to address other species might require investigating environmental responses for those species (e.g., [5,33]).Although it is reasonable to assume, as a first approximation, that the effect of temperature and SWC can be reasonably represented as parabolic functions; it is likely, however, that the relationships are asymmetrical and vary according to tree species (e.g., temperature optimum), soil texture, and available soil water [34].Over periods of 30 days or more, the assumption that photosynthesis (and growth) is a linear function of absorbed PAR is realistic, as long as stomatal conductance is not limiting and leaves are present.However, this assumption may be less realistic at shorter time intervals.iv.The process of normalization, as applied in northwest NB, should be applied to other biomes in a similar manner (i.e., normalization is biome specific).

Comparison of PSD Methods
Features of our method that contribute to the field of PSD modelling (e.g., [12,14,15,18,35]) include: i. Currently, many PSD studies use interpolation to derive input surfaces of biophysical variables in the development of PSD maps.In our work, we employ RS data and process-based modelling techniques to eliminate the inherent problems associated with interpolation of sparse data.ii.The PSD work of Saatchi et al. explicitly uses RS data [18].However, predicting the future state of the landscape under scenarios of climate change would be impossible because their methods lack the capability to predict future states.In contrast, our method can be used to predict future PSD of selected tree species (a) by directly incorporating Global Circulation Model (GCM)-projected increases in temperature in the GDD map, and (b) by replacing estimates of TVWI with SWC derived with the soil-hydrology component of LanDSET.Procedures to carry out such work are addressed in a recent paper [36] discussing the effects of climate change on species distribution in the Acadian forest region of eastern Nova Scotia.iii.Most studies in PSD modelling are applied at spatial resolutions of 1-10 km [e.g., 12,14,15,35].As a result, the utility of the information to forest management and wood supply analysis and planning professionals is limited.Our method can generate realistic maps of PSD at 250-m resolution with existing MODIS data and landscapeprocess modelling techniques.

Concluding Remarks
We have established a simple and practical approach to modelling species-specific PSD for bF by integrating spatial descriptions of biophysical variables derived from RS and process modelling techniques for a largely unmonitored region.The use of RS data to derive estimates of potential species distribution (by way of habitat suitability, HS, calculations) at sufficiently small spatial resolutions (250-m) is a unique contribution to forest management.Direct linkage of environmental variables to plant growth by way of environmental-response functions is an element that permits the current work to be extended to investigate tree distribution and, in the future, growth and yield projections under various scenarios of climate change.The modular form of the definition of HS and PSD makes the approach adaptable to other controlling variables (e.g., soil) and plant species, presenting a foundation for monitoring changes in natural systems.

Figure 1 .
Figure 1.Location of the study area: (a) in northwest New Brunswick, Canada; (b) extentof high bF-content stands (with >50% bF, by volume) from an existing digital forest-cover map (28120 stands, in all).Forest cover in ecoregion 3a is not available for this study.

Figure 2 .
Figure 2. Relation between daily mean a T and T S for the January-October period of 2004 The map was generated by processing 8-day composites of MODIS data of T S at 1-km resolution and surface reflectance at 250-m resolution for the May-October period of 2003-2005.Surface reflectance values were used to generate NDVI for the area at 250-m resolution.Realistic Daily mean air temperature (K

Figure 4 .
Figure 4. Conceptual framework for modelling and validating PSD.
and past management activities (D.E.Swift, personal communication, 2009).Currently, bF resides in the understory and is anticipated to become prominent with the release and growth of bF in the understory as the stand nears maturity (D.E.Edwin Swift, personal communication, 2009).Also, high GDDs in the area potentially favour growth of hardwood species, which would 0