Close Home
iForest - Biogeosciences and Forestry
vol. 9, pp. 178-186
Copyright © 2016 by the Italian Society of Silviculture and Forest Ecology
doi: 10.3832/ifor1712-008

Research Articles

Development of monitoring methods for Hemlock Woolly Adelgid induced tree mortality within a Southern Appalachian landscape with inhibited access

Tuula Kantola (1)Corresponding author, Päivi Lyytikäinen-Saarenmaa (2), Robert N Coulson (1), Markus Holopainen (2), Maria D Tchakerian (1), Douglas A Streett (3)


Eastern hemlock (Tsuga canadensis L. Carr - Pinaceae) and Carolina hemlock (Tsuga caroliniana Engelm. - Pinaceae) are shade tolerant tree species that have a long life cycle ([50], [69]). These species are well adapted to a broad range of soils types and site conditions ([56]). Eastern hemlock has an extensive geographic range in eastern North America and is considered to be a foundation (keystone) species that provides a variety of intermediate and final ecosystem services, e.g., regulation of stream temperature, mitigation of soil erosion ([70]), and habitat and food resources for a suite of taxa, including birds, mammals ([25], [36], [56]), fish, invertebrates, amphibians, and reptiles ([36]). In contrast to eastern hemlock, Carolina hemlock has a more restricted geographic distribution in the Southern Appalachians and is often found on nutrient impoverished soils at higher elevations.

Hemlock woolly adelgid (Adelges tsugae Annand, HWA) is an invasive forest pest that feeds on hemlock parenchyma cells ([74]). The species was introduced from Japan and first detected in the eastern USA in the 1950s ([18]). Initially, the HWA was not considered to be a forest pest in its new environment ([69]). Subsequently, populations of this insect expanded and infested eastern and Carolina hemlock throughout a broad expanse of their ranges ([4]). Herbivory by the HWA results in progressive weakening and eventual mortality of trees ([62]). Hemlocks of all age and size classes are vulnerable ([49]) and infested trees routinely succumb to HWA herbivory, within a span of 5-15 years ([62]).

The effects of HWA herbivory at the landscape-scale ([9]) include widespread elimination of hemlocks. Both species composition and stand structure are fundamentally altered ([14]). One result is a landscape dominated by broad-leaved tree species ([61], [1]). In addition to the landscape-scale effects, responses include changes in transpiration ([10], [15]), carbon cycling ([49], [5]), and energy and nutrient fluxes ([62]).

Hemlocks commonly occur in topographically complex terrain, including steep slopes, deep gorges, ravines, and riparian borders ([47]). Hemlocks are not typically dispersed as large continuous canopies, but rather occur in small patches, often mixed with deciduous or other coniferous tree species ([56], [31]). Furthermore, Young & Morton ([75]) suggested that the patchy nature of eastern hemlock decline may also originate from the influence of landscape-level environmental factors. Landscape patterns and topography can directly influence pest populations and dispersal, or influence indirectly the health and distribution of the host trees ([75]). Traditional ground-based tree- and plot-wise inventory procedures and field surveys for insect-induced foliage alterations and their extent are inadequate. New methods are needed for revealing landscape-scale patterns of hemlock decline within areas that have difficult terrain and sparse road-networks. Remote-sensing methodologies provide an alternative means of assessment that can be both accurate and cost-effective. This approach to assessment provides land-cover maps at various and appropriate spatial and temporal scales ([13], [29]) and, aside from validation, reduces the need for costly and time-consuming ground survey ([43], [45]).

Aerial images are widely used in forest mapping and monitoring, since this medium is commonly available and relatively inexpensive ([26]). Color-infrared images (CIR) with a near-infrared (NIR) band are particularly well suited for tree species recognition, compared with true color images. The NIR band is especially useful for distinguishing conifers from hardwood species ([21]). Low-resolution remote sensing lacks power to reveal fine-scale community structures and dispersion patterns. Furthermore, remote-sensing applications using LiDAR (Light Detection and Ranging) provide another means for accurate mapping of forest vegetation in three dimensions. This approach has a variety of applications in vegetation mapping and monitoring in forestry, and other disciplines ([22]).

Aerially collected images have been widely utilized in forest health-monitoring applications. High-resolution imageries have been successfully applied in monitoring damage by pest insects, including mountain pine beetle (Dendroctonus ponderosae Hopkins - [6], [44], [73]), European spruce bark beetle (Ips typographus L. - [38]), common pine sawfly (Diprion pini L. - [26]), and Jack pine budworm (Choristoneura pinus pinus Free. - [39]). Accompanying advances in detection technologies has been the development of methods for pattern recognition of remotely sensed data ([40], [72]). These methodologies facilitate the use of large volumes of remotely sensed data with high spatial and spectral resolutions.

Despite the acknowledged ecological importance of hemlock species and concern over the decline of hemlock in Southern Appalachian forests, spatially explicit information on the landscape-scale pattern of hemlock mortality is limited. Hemlock is only of marginal economic value in commercial forestry ([4]). The lack of economic interest coupled with the difficulties in field inventory and monitoring tasks may have hindered quantitative evaluation of the ecological impact of hemlock removal in Southern Appalachian forested landscapes.

Methods for monitoring HWA herbivory on the landscape-scale are needed to estimate these impacts. These landscapes are often remote with limited access. Practical approaches include mapping procedures that can be conducted with no or little ground reference. Spatial resolution of the method should be high enough to enable investigation of landscape patterns and topographical relationships, as well as possible pathways for the HWA. The present investigation is a pilot study that aims towards a flexible monitoring system that could also be applied in inaccessible areas. Accordingly, the goal was to employ remote sensing technologies to assess HWA impact on hemlock mortality in Southern Appalachian forest landscapes. The specific objectives include employment of a semi-automatic method for mapping HWA-induced hemlock mortality as well as conifer patches with a living hemlock component within a Southern Appalachian landscape, using high-resolution aerial images and low pulse-density airborne-scanning LiDAR.

Material and methods 

Study area

The approximately 48-km2 wide study area is located in the Grandfather Ranger District (35° 56′ N, 81° 55′ W), Pisgah National Forest, in the Southern Appalachians of western North Carolina (Fig. 1). The landscape vegetation consists primarily of mature, mixed-species forest stands that are characteristic of the ecoregion. The topography is very rugged, with elevations ranging from 1270 m a.s.l. on the ridgetops to 490 m (mean of 880 m) along the Linville River, which flows through the study site. The southwestern part of the area includes an urban component. Approximately half of the study site occurs within the Linville Gorge Wilderness, which has been only modestly logged; less than 5%, according to Newell & Peet ([48]). However, the area has been disturbed, e.g., by chronic acidic deposition and frequent forest fires ([48], [71], [12]). The most recent wildfire occurred in the southwestern part of the study area in October 2000 ([71], [12]). The three major ecological zones that occur within the study site include the Acidic Cove, Xeric Pine-Oak Heath and Oak Heath, and Mesic Oak-Hickory zones ([60], [1]). This structurally diverse landscape and humid climate provides a variety of habitats for more than 400 vascular plants and a rich diversity of tree species ([59], [53], [60]). A sparse road-network resulting from the complex terrain and preserved wilderness area further inhibits access and data collection in the area. Both eastern and Carolina hemlock are abundant within the area ([24]). The HWA was first detected in the area in the early 2000s ([32]) and was the most significant tree mortality agent during the time period of the study. Patches of dead hemlocks can be observed in the landscape. Hemlock decline is the most abundant on the riverside and along the drainages. We also found that the suppressed small hemlocks growing on the Linville River hillsides were infested and heavily defoliated in May 2014 (Kantola & Lyytikäinen-Saarenmaa, field observation, May 9-11, 2014).

Fig. 1 - The study area. Location of the Grand Father Ranger District (green) and the study area (red) in the western North Carolina (©ESRI - left) and an aerial image mosaic of the study area (right).

Remote sensing data sets

Aerial mages

CIR imagery with 1-m spatial resolution was acquired with a multiple Intergraph Digital Mapping Camera system (Intergraph Corp., Huntsville, AL, USA) at an altitude of approximately 9000 m in the summer of 2012 (leaf-on) by the National Agricultural Inventory Program (NAIP). The images were captured simultaneously from four 3072 × 2048 pixel multispectral cameras with 30 mm lenses producing red, green, blue, and NIR bands. The CIR imagery used in the analysis was resized to 3 m × 3 m pixel size to better correspond to the size of an individual tree crown.

Color aerial imagery (RGB) with 15-cm spatial resolution in the red, green, and blue bands was acquired by the Sanborn Map Company Inc. (Colorado Springs, CO, USA) with a large-format Zeiss (Carl Zeiss AG, Oberkochen, Germany) / Intergraph DMC in the winter 2010 (leaf-off). The flying altitude was approximately 1500 m above the mean terrain. The RGB imagery was used as a reference to enhance the accuracy of the training and testing data sets.

Airborne-scanning LiDAR

Airborne-scanning LiDAR was acquired during the North Carolina Floodplain Mapping Program phase II, in 2003. The LiDAR point cloud was acquired with a Leica GeoSystems Aeroscan system (Leica Geosystem AG, Heerbrugg, Switzerland). The flying altitude was 3000 m above the mean terrain at a speed of 150 knots, with a field of view of 55 degrees, and a laser pulse rate of 23100 Hz. The density of the pulses returned within the area was less than 1 per m2. A LiDAR point cloud of approximately 40 km2 was downloaded from the U.S. Geodetic Survey (USGS) Earth Explorer ([66]). The area covered by the LiDAR point cloud was used in the analysis. A canopy height model (CHM) was derived from the LiDAR point cloud and used in the first sub-task of the study.


General workflow

The first sub-task was to exclude non-forested land from the study site, including urban areas, water, and bare ground. We also excluded deep black shadows that lacked spectral information. These tasks were accomplished by creating a forest mask, using CHM and Normalized Difference Vegetation Index (NDVI) layers. The rationale for the forest mask creation was to avoid overestimation of dead tree coverage due to similar spectral reflectance associated with some non-forested patch types, such as bare ground and roads. Furthermore, reducing the number of cover classes usually improves the overall classification accuracy. The second sub-task included the classification of the remaining forest canopy cover. We assessed the accuracy by calculating the classification accuracies and Cohen’s kappa-values (κ - [7]) for both sub-tasks separately. The implementation was done by employing the FUSION® software (FUSION/LDV, USDA Forest Service, Seattle, WA, USA - [42]), Environment for Visualizing Images (Exelis Visual Information Solutions - ENVI®) software (EXELISVIS Inc., Boulder, CO, USA), and the ArcGIS® Geographic Information System environment (Environmental Systems Research Institute - ESRI, Redlands, CA, USA). The study workflow was divided into four main steps as follows:

  1. creation of training and testing data sets via visual interpretation (section “Training and testing data sets”);
  2. creation of LiDAR-derived CHM and NDVI layers (section “Forest mask creation”);
  3. forest mask creation with decision-tree classification and extraction of non-forested areas from the aerial imagery (section “Forest mask creation”); and
  4. conducting Support Vector Machine (SVM) classification for the forest cover classes (section “Forest cover classification”).

Training and testing data sets

We created testing and training data sets by visual interpretation of the aerial images. The original images (CIR and RGB) were used in the visual interpretation. The corresponding pixels from the resized image mosaic were chosen for the training and testing data sets. A 200 × 200-m systematic point network was created to test the accuracy of the forest mask created in sub-task 1. A total of 1080 pixels were assigned to the forest and non-forest classes, yielding 792 forest and 288 non-forest pixels. We created separate training and testing data sets for the second sub-task to classify and evaluate the forest cover classes (conifers, hardwood species, and dead trees). A sample of pixels was chosen that could be classified by an expert, without high uncertainty, into conifers, hardwood species, and dead trees. A total of 7925 pixels (5701 for training and 2224 for testing) were assigned to these data sets. A total of 81 041 m2 of the study area were used for training and testing (9720 m2 for the forest mask, 71 321 m2 for the forest cover classification - Tab. 1). The testing data set covered approximately 40% of the reference data set.

Tab. 1 - Areas and number of pixels in training and testing data sets for three different forest cover classes.

Forest mask creation

Separation of above-ground surface features offers useful data for a wide range of environmental applications. The digital surface model is a three-dimensional representation of the ground surface in non-vegetated terrain, as well as aboveground features, such as vegetation and buildings. In forestry applications, they are usually referred to as CHMs. The CHM was created from the LiDAR point cloud via the following steps. The elevation of the highest return within each grid cell was assigned as the local maximum to the grid cell center. The ground elevation was estimated, using a ground filter algorithm described in detail by Kraus & Pfeifer ([35]). Laser heights above ground (normalized heights, or canopy heights) were then calculated by subtracting the ground elevations from the corresponding local maxima. The CHM layer created with a spatial resolution of 3 × 3-m was smoothened with a 3 × 3 pixel median filter.

The NDVI is widely used in ecological applications ([54]). The principle behind it is that chlorophyll absorbs visible light in the read region of the electromagnetic spectrum, whereas radiation in the nonvisible NIR region of the spectrum is scattered by the leaf structure of a plant ([46]). As a result, vigorously growing healthy vegetation has low red-light reflectance and high NIR reflectance and, hence, high NDVI values.

Decision-tree classification is a non-parametric procedure and represents a practical method for land cover classification ([64]). A decision tree can be defined as a classification procedure that repetitively partitions a data set into smaller subdivisions ([16]). The partition is made, based on tests defined at each node of the decision tree. Decision-tree classification has advantages over traditional supervised classification methods, such as maximum likelihood or neural networks, which are widely used in remote sensing applications ([40], [64]). There are no required assumptions about the distributions and variations in data. The procedure accommodates missing values and the use of both numerical and categorical input values ([16], [40]). The method is flexible and can also process non-linear relationships ([40]).

We used the CHM and the NDVI layers in the decision tree classification for delineation of the forested areas. CHM heights greater than 1.2 m were considered vegetation, buildings, and other objects above the ground. A typical threshold height for vegetation returns in forestry application is 2 m ([26]). We kept the threshold value for the canopy height even lower, because the LiDAR point cloud was somewhat dated. With a 1.2-m threshold, grasslands and other short vegetation could still be excluded. This concession was reasonable, because we were simply separating terrain and vegetation and not estimating tree or stand characteristics. A low threshold value of 0.05 for the NDVI was chosen. We excluded deep shadows, water, and other non-vegetation elements having low NDVI values as well. The NDVI threshold chosen, however, did not exclude the dead trees. The resulting classification layer was filtered with a neighborhood majority filter to smooth the forest mask. The area within the forest mask was further used in the second sub-task.

Forest cover classification

The goal of the forest cover classification was to detect patches of dead hemlocks and distinguish living conifers from hardwood species. The procedure permits mapping of clusters of dead trees and potential areas of living hemlocks. The spectral reflectance of deciduous tree species is known to differ from that of conifers. The SVM classification methods presented by Boser et al. ([3]) and Vapnik ([67]), are not as well-known as many other procedures ([45]). However, their performance can exceed that of other classification methods ([17], [45]). SVM is a supervised non-parametric learning technique ([8]). SVM aims to determine the location of optimal decision boundaries separating different classes ([67]). The nearest data points to the resulting hyperplane that are used to measure the margin are called support vectors ([52]). SVM approaches using kernel functions can map non-linear data into a higher dimensional space. In that space, a linear separating surface between two classes is searched ([17]). SVM classifiers are suitable for remote sensing applications, especially with limited training data sets ([41], [45]). SVM classifiers are also seen as being robust to noise in data sets and to high-dimensional data. SVM minimizes classification errors without any prior assumptions about the probability distributions of a data set.

SVM classifiers have been successfully used in remote sensing applications, including vegetation classification ([17], [37]), forest classification ([23]), and tree species classification ([19]). The radial basis function (RBF) and polynomial kernels are the most commonly used in remote sensing-based SVM classifications ([29]). Kavzoglu & Colkesen ([29]) gained better overall classification accuracy with RBF in land cover classification. We used the spectral bands of red, green, NIR, and NDVI layers, and the RBF kernel in the SVM classification.


Forest mask

We extracted the non-forested areas from the CIR image for the second phase of the analysis (Fig. 2). The area of the LiDAR point cloud was smaller than the CIR mosaic causing the black frame. The overall classification accuracy was 93.5% (κ = 0.84).

Fig. 2 - The used high resolution image tiles. The original color-infrared (CIR) image mosaic of the study area (left), and the classified image for the second phase (right). On the classified image, only the forested areas are visible and the non-forested areas are extracted. Black color indicates the extracted areas.

Forest cover classification

The forested area, a total of 30.2 km2, was classified into three forest cover classes: dead trees, conifers, and hardwood. The overall classification accuracy was 98.1% (κ = 0.96). Hardwood species covered over 55% of the classified area, while conifers and dead trees covered 44.7% of the classified area (42.6% and 2.1%, respectively) (Tab. 2). Over 0.6 km2 of the area was classified as dead trees and 12.9 km2 as conifers. The proportions found throughout the study area, including unclassified areas, were 1.63% dead trees, 42.5% hardwood species, 32.8% conifers, and 23.1% non-forested/non-classified (Fig. 3).

Tab. 2 - (A) Classification accuracies within the study area: an error confusion matrix of the SVM classification results. (B) The area covered by different forest cover classes and the proportions of the forested area.
Fig. 3 - The resulting classification image. The SVM classification image (left): 1=dead trees; 2=hardwood species; 3=conifers; and 4= unclassified. A magnified portion of the classification image (a black square) is illustrated on the left side of the figure. The classification image is on top and the CIR image below (right). Conifers show as darker purple color in the CIR image, hardwood species show as pink or red, and dead trees show as pale green.

The classification image combined with the digital terrain model (DTM) revealed that conifers are most abundant in drainages and on northern and western slopes (Fig. 4). Clusters of dead trees could be found, especially near the Linville River, in drainages and at higher elevations.

Fig. 4 - The classification result combined with a digital terrain model at the Linville River Gorge. The classification image pronounces areas under interest, conifers and dead trees (green and red, respectively). Pixels classified as non-forest or hardwood species show as grey.


Evaluating hemlock mortality

Previous studies directed at evaluating the impacts of hemlock decline have been based mainly on plot-wise field assessments and the use of lower-resolution remote sensing approaches. Plot-wise sampling schemes have been used for studying both small-scale within-stand effects of hemlock mortality ([50], [11], [34]) and broad-scale regional and state-wide impacts ([65]). Estimates for rates of hemlock mortality varied considerably. Results of plot- and stand-wise studies revealed hemlock mortality rates between 0% and 95%, depending on the infestation histories and latitudes ([50], [11], [34]).

Landsat Thematic Mapper (TM) has been the most commonly used remote sensing sensor in the mapping of hemlocks and HWA-induced hemlock mortality ([2], [58], [71], [33]). These studies provided good insights into HWA herbivory and hemlock decline at different spatial scales. On the other hand, wall-to-wall studies at landscape level with high spatial resolution are scarce. Kantola et al. ([28]) studied spatial pattern of hemlock decline in the Southern Appalachians (NC, USA). They visually assessed tree mortality in the upper canopy cover layer within the Lower Linville Watershed. They suggested that despite the modest loss in total biomass, corresponding to 0.1% of the canopy surface, the impacts were substantial to the area due to the patchy nature of the hemlock decline, as well as the elimination of a foundation species. However, on the landscape scale, the magnitude of impacts is still mainly unknown and the lack of information is even more pronounced in inaccessible areas.

Hemlocks are distributed throughout forested landscapes and grow in mixed-species stands. These landscapes generally comprise several tree species, including other conifers. High-resolution remote sensing enables a more detailed mapping of broad spatial extent. The high-resolution imagery used in our study enable a more detailed projection of the area and pattern of the hemlock decline suitable for estimating the magnitude of ecological and social impacts on the landscape-scale ([9]).

Two-phase forest cover classification

The first phase of the classification comprised the creation of a forest mask. This sub-task was mandatory for distinguishing dead trees from bare ground and some urban infrastructure, such as roads. Without this phase, our assessment for dead trees would have been overestimated. We observed that the mask followed the borders of forest and non-forested in most parts. An overall classification accuracy of 93.5% (κ = 0.83) was obtained for forest/ non-forest classification. Minor inaccuracies could be found, e.g., on the edges of bare ground ridges. This error may have been partly due to filtering of the mask and loss of small-scale information.

Projection transformations and photogrammetry can also affect spatial accuracy and thereby result in imprecise overlapping of the remote sensing data sets. The LiDAR point cloud used in the study was acquired in 2003, nine years earlier than the aerial image mosaic. A destructive forest fire occurred in the study area in 2000 ([12]), which certainly affected the forest stand dynamics. There is no information available on logging outside of the Linville Wilderness area. The threshold value for the height was purposefully kept low to minimize the impact of changes in the vegetation between 2003 and 2012. More recent LiDAR data could have improved the accuracy of the method. With simultaneously acquired high-pulse density LiDAR and aerial imagery, it may be possible to delineate and classify individual trees ([26]).

The results of our forest cover classification indicated that dead trees occurred in an area of approximately 0.64 km2 (about 2% of the study area). We found that 12.9 km2 of the area were covered by conifers. The last high mortality-causing disturbance occurred in the early 2000s. The area was infested with the southern pine beetle (Dendroctonus frontalis Zimm., SPB - [30]). Inspection of a time-series of high resolution aerial images, acquired from the Linville Gorge revealed that a large portion of the dead trees in the forest canopy surface became invisible within a span of five or more years after dying (Kantola et al., unpublished data). We assume that the majority of the dead trees were hemlocks, because the HWA has been the only high mortality agent during recent years, Therefore, we suggest that HWA herbivory has been the cause of death to nearly 5% of the overstory coniferous component of the vegetation community present at the time of data acquisition. Other plausible tree mortality agents within the area include minor infestations of other species and abiotic factors, including hardwood species. Kantola et al. ([28]) found that dead trees covered an area approximately of 0.72 km2 (0.1%) of the forested area within the Lower Linville Watershed. The central parts of the area in our study belong to the same watershed.

Use of passive remote sensing, i.e., aerial images enables only investigation of the overstory HWA mortality. Orwig & Foster ([50]) and Krapfl et al. ([34]) suggest that the mortality may be greater among understory hemlocks than in the upper canopy layers. However, the information provided by aerial images can produce good estimates of the extent and pattern of hemlock decline that include most of the biomass.

We obtained an overall accuracy of 98.1% (κ = 0.96) for the three forest cover classes. These values are probably an overestimation resulting from the visual assessment of aerial images. In cases where there was uncertainty, pixels were excluded from the evaluation data set. Shadows were excluded from the forest cover classification. Topographic correction methods for high-resolution images in order to reduce the effect of shadows can be problematic. With planning, the image acquisition can be done on date and at time of day when the shadows are less pronounced. Another option is to acquire two sets of images captured at different times of day for the assessment. Most of the shadowed areas were assumed to be covered by forest. Excluding the areas may have introduced bias into the results. However, our classification approach provided reasonable accuracy for mapping fine-scale HWA-induced tree mortality. Thomlinson et al. ([63]) set the criteria for successful land-cover classifications stating that overall classification accuracy should exceed 85% with at least 70% per-class accuracy.

A large component of the coniferous community within the study area was assumed to be hemlock. Several sources support this assumption. Using the North Carolina Vegetation Survey (NCVS) protocol, Newell & Peet ([48]) sampled 181 plots within the Linville River Watershed. Eastern and Carolina hemlocks were abundant in the vegetation classes covering most of the northern part of the Linville Gorge Wilderness, which included much of the study area ([48]). Knebel & Wentworth ([30]) reported that a combination of frequent forest fires and previous SPB infestations diminished the pine (Pinus spp.) component of the vegetation community. In addition to eastern and Carolina hemlocks, possible coniferous species within the study area include eastern white pine (P. strobus L.), pitch pine (P. rigida P.Mill.), table mountain pine (P. pungens Lamb.), and Virginia pine (P. virginiana Mill. - [48], [12]).

Since hemlocks and other conifers have similar reflectance values, distinguishing among the species is difficult and prone to errors ([57], [51], [31]). One option for addressing this issue is to sample accessible parts of target areas that are classified as conifers and collect tree-wise field data with accurate locations. Reference data may be collected from another, and accessible, landscape with similar tree species composition, with consideration. A spectral library for the coniferous species present could be built using similar remote sensing data. These spectral reflectance signature may help to distinguish the hemlock component from other conifers. Use of existing auxiliary information under conditions typical of hemlock sites could be another approach to separate hemlocks and other conifers. Environmental layers, such as topography, soil, temperature, and site-type layers can be utilized. A vast array of algorithms could be employed in modeling of the most probable areas for hemlock species within the conifer patches.

An early symptom of HWA infestation is gradual defoliation ([50], [11]). Defoliation can be estimated to some degree from aerial images or LiDAR point clouds ([55], [27], [68]). For example, Elliott & Vose ([11]) discovered that after less than three years of an HWA infestation, the mean defoliation level was over 80% in the sampling plots with 100% hemlock infestation rate. In the present study, defoliated hemlocks were included in the forest cover class of conifers. A portion of the heavily defoliated trees may have been classified as dead trees.

Remote sensing in hemlock mapping and HWA monitoring

Although high-resolution imageries have not been utilized in mapping hemlocks and HWA-induced tree decline, there are studies that have used various remote sensing approaches. Koch et al. ([31]) studied hemlock mapping via ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer) satellite images in the Great Smoky Mountains of North Carolina. Auxiliary raster layers were used to enhance the classification. Their overall classification accuracy was 85.3% (κ=0.77) and the accuracy of hemlock stand detection was 69%. In the western part of their study area, hemlocks were limited to narrow riparian corridors. In the eastern part, hemlocks were more broadly distributed. These authors assumed that restricting hemlocks to riparian corridors in the training data caused misclassification in the other areas.

Royle & Lathrop ([58]) investigated the use of multi-temporal Landsat TM derived vegetation indexes in mapping the health of hemlock forest stands in Connecticut. The best overall classification accuracy for four hemlock forest health classes was 82%. Wimberly & Reilly ([71]) used Landsat TM images to assess fire severity and species diversity in the Linville River Gorge. Eastern hemlock was associated with moist areas along ravines and at the bottom of the gorge. Kong et al. ([33]) used Landsat TM and ASTER images combined with environmental variables to detect hemlock stands in a study conducted in Kentucky. From the evergreen pixels, they correctly detected 72% of the hemlock points. Pontius et al. ([55]) detected hemlock-dominated pixels from AVIRIS hyperspectral images (Airborne Visible/Infrared Imaging Spectrometer) with 83% accuracy in New York.

Impacts of HWA on the landscapes

Our classification image indicated that conifers and dead trees were abundant in the proximity of the Linville River, on steep hillsides, and at high elevations. Most of the dead trees were located in the northern and northwestern aspects. Kantola et al. ([28]) obtained similar results, showing that within the Lower Linville Watershed, the density of the dead trees was higher in proximity to the Linville River, at higher elevations, and in the northern and northwestern aspects. They also found that the spatial pattern of the dead trees was typically clustered. They obtained the results with more time-consuming visual assessment of aerial images. The classification image showed that conifers were most abundant on northern and western slopes and drainage. This observation is in line with the results from other reports. Narayanaraj et al. ([47]) examined the topographical pattern of living eastern hemlock in the Coweeta Basin of North Carolina. Their results showed that living eastern hemlocks were more abundant close to streams and on flat or gentle slopes, at lower elevations. They found that eastern hemlock was absent above 1250 m. In a study conducted in central Connecticut, Orwig et al. ([51]) observed that a large number of hemlock stands were located on ridge tops, steep hillsides, and in narrow valleys. Hemlock mortality occurred mostly on the westward-facing slopes. Several studies suggested that microclimate and soil conditions related to variation in topography are important factors in HWA damage ([20] [47]).

Trotter et al. ([65]) used the USDA Forest Service Forest Inventory and Analysis (FIA) database to address changes in hemlock abundance in the eastern USA. They observed an increase in the basal area of living and dead hemlocks during a 20-year study period at county- and state-levels. Trotter et al. ([65]) suggested that impacts of HWA are not evident at the regional level. Our results indicate that the reduction in total plant biomass within the landscape is modest. Hemlocks, serving as foundation species, play multiple functional roles in forested landscapes, such as modifying environmental conditions and providing forage and habitat resources for a variety of taxa. Therefore, the impacts of the decline are profound. The clustered pattern of the dead trees, especially in the riparian areas, can intensify the impacts on forest landscapes and ecosystems.

The developed two-phase classification strategies used in this investigation can be adapted for monitoring other invasive pest insects that cause heavy defoliation and tree mortality, e.g., the gypsy moth (Lymantria dispar L.), SPB, and emerald ash borer (Agrilus planipennis Fairmaire). The methodology can also be used with auxiliary information to produce training sets for large area inventories with lower-resolution remote sensing data sets.


The HWA is a major mortality agent of eastern and Carolina hemlocks throughout their distribution in the eastern USA. Accurate spatially explicit information on hemlock distribution and HWA-induced mortality is not available at the landscape-level. Therefore, it is not possible to evaluate the impact of this insect on the vegetation community in which hemlocks are a prominent component. The straightforward two-phase methodology described herein, which utilizes aerial images and LiDAR point clouds, can result in distribution maps of considerable accuracy. This procedure is suitable for mapping the coniferous patches that include the hemlock component of the vegetation community and for distinguishing tree mortality resulting from the HWA or other disturbance agents.


This study was made possible by financial aid from the Graduate School in Forest Sciences (GSForest) in Finland, the Finnish Academy project “Centre of Excellence in Laser Scanning Research” (CoE-LaSR, decision number 272195), the University of Helsinki project “Towards semi-supervised characterization and large-area planning of forest resources using airborne laser scanning data acquired for digital elevation modeling”, and by the US Forest Service through USDA Forest Service cooperative agreement SRS-12-CA-11330129-077. We thank Dr. Hannu Saarenmaa for fruitful insights.


Birt AG, Zeng Y, Tchakerian MD, Coulson RN, Lafon DM, Cairns DM, Waldron J, Xi W, Chen S-H, Street DA (2014). Evaluating Southern Appalachian forest dynamics without eastern hemlock: consequence of herbivory by the hemlock woolly adelgid. Open Journal of Forestry 4 (2): 91-98.
::CrossRef::Google Scholar::
Bonneau LR, Shields KS, Civco DL (1999). Using satellite images to classify and analyze the health of hemlock forests infested by the hemlock woolly adelgid. Biological Invasions 1 (2-3): 255-267.
::CrossRef::Google Scholar::
Boser BE, Guyon IM, Vapnik VN (1992). A training algorithm for optimal margin classifiers. In: Proceedings of the “5th annual workshop on computational learning theory” (Haussler D ed). Pittsburgh (PA, USA) 27-29 July 1992. ACM Press, Pittsburgh, PA, USA, pp. 144-152.
::CrossRef::Google Scholar::
Clark JT, Fei S, Liang L, Rieske LK (2012). Mapping eastern hemlock: comparing classification techniques to evaluate susceptibility of a fragmented and valued resource to an exotic invader, the hemlock woolly adelgid. Forest Ecology and Management 266: 216-222.
::CrossRef::Google Scholar::
Cobb RC (2010). Species shift drives decomposition rates following invasion by hemlock woolly adelgid. Oikos 119 (8): 1291-1298.
::CrossRef::Google Scholar::
Coggins S, Coops NC, Wulder MA (2008). Initialization of an insect infestation spread model using tree structure and spatial characteristics derived from high spatial resolution digital aerial imagery. Canadian Journal of Remote Sensing 34 (6): 485-502.
::CrossRef::Google Scholar::
Cohen JA (1960). Coefficient of agreement for nominal scales. Educational and Psychological Measurement 20: 37-46.
::CrossRef::Google Scholar::
Cortes C, Vapnik V (1995). Support-vector networks. Machine Learning 20 (3): 273-297.
::Google Scholar::
Coulson RN, Tchakerian MD (2010). Basic landscape ecology. KEL Partners Incorporated, Boston, USA, pp. 300.
::Online::Google Scholar::
Daley MJ, Phillips NG, Pettijohn C, Hadley JL (2007). Water use by eastern hemlock (Tsuga canadensis) and black birch (Betula lenta): implications of effects of the hemlock woolly adelgid. Canadian Journal of Forest Research 37 (10): 2031-2040.
::CrossRef::Google Scholar::
Elliott KJ, Vose JM (2011). The contribution of the Coweeta Hydrologic Laboratory to developing an understanding of long-term (1934-2008) changes in managed and unmanaged forests. Forest Ecology and Management 261: 900-910.
::CrossRef::Google Scholar::
Elliot KJ, Knoepp JD, Vose JM, Jackson WA (2013). Interacting effects of wildfire severity and liming on nutrient cycling in a Southern Appalachian wilderness area. Plant and Soil 366 (1-2): 165-183.
::CrossRef::Google Scholar::
Foody GM (2002). Status of land cover classification accuracy assessment. Remote Sensing of Environment 80 (1): 185-201.
::CrossRef::Google Scholar::
Ford CR, Elliott KJ, Clinton BD, Kloeppel BD, Vose JM (2012). Forest dynamics following eastern hemlock mortality in the Southern Appalachians. Oikos 121 (4): 523-536.
::CrossRef::Google Scholar::
Ford CR, Vose JM (2007). Tsuga canadensis (L.) Carr. mortality will impact hydrologic processes in Southern Appalachian forest ecosystems. Ecological Applications 17 (4): 1156-1167.
::CrossRef::Google Scholar::
Friedl MA, Brodley CE (1997). Decision tree classification of land cover from remotely sensed data. Remote Sensing of Environment 61 (3): 399-409.
::CrossRef::Google Scholar::
Gualtieri JA, Cromp RF (1999). Support vector machines for hyperspectral remote sensing classification. In: Proceedings of the “27th AIPR Workshop: Advances in Computer Assisted Recognition” (Merisko RJ ed). Washington (DC, USA) 27 Oct 1999. SPIE, Washington, DC, USA, pp. 221-232.
::CrossRef::Google Scholar::
Havill NP, Montgomery ME, Yu G, Shiyake S, Caccone A (2006). Mitochondrial DNA from hemlock woolly adelgid (Hemiptera: Adelgidae) suggests cryptic speciation and pinpoints the source of the introduction to eastern North America. Annals of the Entomological Society of America 99 (2): 195-203.
::CrossRef::Google Scholar::
Heikkinen V, Tokola T, Parkkinen J, Korpela I, Jääskelainen T (2010). Simulated multispectral imagery for tree species classification using support vector machines. IEEE Transactions on Geoscience and Remote Sensing 48 (3): 1355-1364.
::CrossRef::Google Scholar::
Hodkinson ID (2005). Terrestrial insects along elevation gradients: species and community responses to altitude. Biological Reviews 80 (3): 489-513.
::CrossRef::Google Scholar::
Holmgren J, Persson A (2004). Identifying species of individual trees using airborne laser scanner. Remote Sensing of Environment 90 (4): 415-423.
::CrossRef::Google Scholar::
Holopainen M, Vastaranta M, Liang X, Hyyppä J, Jaakkola A, Kankare V (2014). Estimation of forest stock and yield using Lidar data. In: “Remote Sensing of Natural Resources” (Wang G, Weng Q eds). CRC Press, Taylor and Francis Group, Boca Raton, FL, USA, pp. 259-290.
::Online::Google Scholar::
Huang H, Gong P, Clinton N, Hui F (2008). Reduction of atmospheric and topographic effect on Landsat TM data for forest classification. International Journal of Remote Sensing 29 (19): 5623-5642.
::CrossRef::Google Scholar::
Jetton RM, Dvorak WS, Whittier WA (2008). Ecological and genetic factors that define the natural distribution of Carolina hemlock in the southeastern United States and their role in ex situ conservation. Forest Ecology and Management 255: 3212-3221.
::CrossRef::Google Scholar::
Jordan JS, Sharp WM (1967). Seeding and planting hemlock for ruffed grouse cover. Research paper NE-83, USDA Forest Service, Upper Darby, PA, USA, pp. 17.
::Online::Google Scholar::
Kantola T, Vastaranta M, Yu X, Lyytikäinen-Saarenmaa P, Holopainen M, Talvitie M, Kaasalainen S, Solberg S, Hyyppä J (2010). Classification of defoliated trees using tree-level airborne laser scanning data combined with aerial images. Remote Sensing 2: 2665-2679.
::CrossRef::Google Scholar::
Kantola T, Vastaranta M, Lyytikäinen-Saarenmaa P, Holopainen M, Kankare V, Talvitie M, Hyyppä J (2013). Classification of needle loss of individual Scots pine trees by means of airborne laser scanning. Forests 4 (2): 386-403.
::CrossRef::Google Scholar::
Kantola T, Lyytikäinen-Saarenmaa P, Coulson RN, Strauch S, Tchakerian MD, Holopainen M, Saareenmaa H, Streett DA (2014). Spatial distribution of hemlock woolly adelgid induced hemlock mortality in the Southern Appalachians. Open Journal of Forestry 4 (05): 492-506.
::CrossRef::Google Scholar::
Kavzoglu T, Colkesen I (2009). A kernel functions analysis for support vector machines for land cover classification. International Journal of Applied Earth Observation and Geoinformation 11 (5): 352-359.
::CrossRef::Google Scholar::
Knebel L, Wentworth TR (2007). Influence of fire and southern pine beetle on pine-dominated forests in the Linville Gorge Wilderness, North Carolina. Castanea 72 (4): 214-225.
::CrossRef::Google Scholar::
Koch FH, Cheshire HM, Devine HA (2005). Mapping hemlocks via tree-based classification of satellite imagery and environmental data. Forest Health Technology Enterprise Team 2005-01, USDA Forest Service, Asheville, NC, USA, pp. 109-115. -
::Online::Google Scholar::
Koch FH, Chesire HM, Devine HA (2006). Landscape-scale prediction of hemlock woolly adelgid, Adelges tsugae (Homoptera: Adelgidae), infestation in the southern Appalachian Mountains. Environmental Entomology 35 (5): 1313-1323.
::CrossRef::Google Scholar::
Kong N, Fei S, Rieske-Kinney L, Obrichy J (2008). Mapping hemlock forest in Harlan County, Kentucky. In: Proceedings of the “6th Southern Forestry and Natural Resources GIS Conference” (Bettinger P, Merry P, Fei K, Drake S, Nibbelink J, Hepinstall N, Athens J eds). Orlando (FL, USA) 24-26 Mar 2008. Warnell School of Forestry and Natural Resources, University of Georgia, Athens, GA, USA, pp. 107-117.
::Online::Google Scholar::
Krapfl KJ, Holzmueller EJ, Jenkins MA (2011). Early impacts of hemlock woolly adelgid in Tsuga canadensis forest communities of the southern Appalachian Mountains. Journal of the Torrey Botanical Society 138: 93-106.
::CrossRef::Google Scholar::
Kraus K, Pfeifer N (1998). Determination of terrain models in wooded areas with airborne laser scanner data. ISPRS Journal of Photogrammetry and Remote Sensing 53 (4): 193-203.
::CrossRef::Google Scholar::
Lapin B (1994). The impact of hemlock woolly adelgid on resources in the Lower Connecticut River Valley. Report for the NE Center for Forest Health Research, USDA Forest Service, Hamden, CT, USA, pp. 43.
::Google Scholar::
Lardeux C, Frison PL, Tison C, Souyris JC, Stoll B, Fruneau B, Rudant JP (2009). Support vector machine for multifrequency SAR polarimetric data classification. IEEE Transactions on Geoscience and Remote Sensing 47 (12): 4143-4152.
::CrossRef::Google Scholar::
Lausch A, Heurich M, Gordalla D, Dobner HJ, Gwillym-Margianto S, Salbach C (2013). Forecasting potential bark beetle outbreaks based on spruce forest vitality using hyperspectral remote-sensing techniques at different scales. Forest Ecology and Management 308: 76-89.
::CrossRef::Google Scholar::
Leckie DG, Cloney E, Joyce SP (2005). Automated detection and mapping of crown discolouration caused by jack pine budworm with 2.5 m resolution multispectral imagery. International Journal of Applied Earth Observation and Geoinformation 7 (1): 61-77.
::CrossRef::Google Scholar::
Mahesh P, Mather PM (2003). An assessment of the effectiveness of decision tree methods for land cover classification. Remote Sensing of Environment 86 (4): 554-565.
::CrossRef::Google Scholar::
Mantero P, Moser G, Serpico SB (2005). Partially supervised classification of remote sensing images through SVM-based probability density estimation. IEEE Transactions on Geoscience and Remote Sensing 43 (3): 559-570.
::CrossRef::Google Scholar::
McGaughey RJ (2009). FUSION/LDV: software for LiDAR data analysis and visualization. USDA Forest Service, Pacific Northwest Research Station, Seattle, WA, USA, pp. 123.
::Google Scholar::
Means JE, Acker SA, Fitt BJ, Renslow M, Emerson L, Hendrix CJ (2000). Predicting forest stand characteristics with airborne scanning LiDAR. Photogrammetric Engineering and Remote Sensing 66 (11): 1367-1372.
::Online::Google Scholar::
Meddens AJ, Hicke JA, Vierling LA (2011). Evaluating the potential of multispectral imagery to map multiple stages of tree mortality. Remote Sensing of Environment 115 (7): 1632-1642.
::CrossRef::Google Scholar::
Mountrakis G, Im J, Ogole C (2011). Support vector machines in remote sensing: a review. ISPRS Journal of Photogrammetry and Remote Sensing 66 (3): 247-259.
::CrossRef::Google Scholar::
Myneni RB, Hall FG, Sellers PJ, Marshak AL (1995). The interpretation of spectral vegetation indexes. IEEE Transactions on Geoscience and Remote Sensing 33 (2): 481-486.
::CrossRef::Google Scholar::
Narayanaraj G, Bolstad PV, Elliott KJ, Vose JM (2010). Terrain and landform influence on Tsuga canadensis (L.) Carrière (eastern hemlock) distribution in the southern Appalachian Mountains. Castanea 75: 1-18.
::CrossRef::Google Scholar::
Newell CL, Peet RK (1998). Vegetation of Linville Gorge Wilderness, North Carolina. Castanea 63 (3): 275-322.
::Online::Google Scholar::
Nuckolls AE, Wurzenburger N, Ford CR, Hendrick RL, Vose JM, Kloeppel BD (2009). Hemlock declines rapidly with hemlock woolly adelgid infestation: impacts on the carbon cycle of Southern Appalachian forests. Ecosystems 12: 179-190.
::CrossRef::Google Scholar::
Orwig DA, Foster DR (1998). Forest response to the introduced hemlock woolly adelgid in southern New England, USA. Journal of the Torrey Botanical Society 125: 60-73.
::CrossRef::Google Scholar::
Orwig DA, Foster DR, Mausel DL (2002). Landscape patterns of hemlock decline in New England due to the introduced hemlock woolly adelgid. Journal of Biogeography 29 (10- 11): 1475-1487.
::CrossRef::Google Scholar::
Pal M, Mather PM (2005). Support vector machines for classification in remote sensing. International Journal of Remote Sensing 26 (5): 1007-1011.
::Online::Google Scholar::
Peet RK, Wentworth TR, White PS (1998). A flexible, multipurpose method for recording vegetation composition and structure. Castanea 63 (3): 262-274.
::Online::Google Scholar::
Pettorelli N, Vik JO, Mysterud A, Gaillard JM, Tucker CJ, Stenseth NC (2005). Using the satellite-derived NDVI to assess ecological responses to environmental change. Trends in Ecology and Evolution 20 (9): 503-510.
::CrossRef::Google Scholar::
Pontius J, Hallett R, Martin M (2005). Using AVIRIS to assess hemlock abundance and early decline in the Catskills, New York. Remote Sensing of Environment 97 (2): 163-173.
::CrossRef::Google Scholar::
Quimby JW (1995). Value and importance of hemlock ecosystems in the eastern United States. In: Proceedings of the “First Hemlock Woolly Adelgid Review” (Salom S, Tigner T, Reardon RC eds). Forest Health Technology Enterprise Team 96-10, USDA Forest Service, Morgantown, WV, USA, pp. 1-8.
::Online::Google Scholar::
Royle DD, Lathrop RG (1997). Monitoring hemlock forest health in New Jersey using Landsat TM data and change detection techniques. Forest Science 43 (3): 327-335.
::Online::Google Scholar::
Royle DD, Lathrop RG (2002). Discriminating Tsuga canadensis hemlock forest defoliation using remotely sensed change detection. Journal of Nematology 34 (3): 213-221.
::Online::Google Scholar::
Schafale MP, Weakly AS (1990). Classification of the natural communities of North Carolina: third approximation. NC Natural Heritage program, Division of Parks and Recreation, Raleigh, NC, USA, pp. 321.
::Online::Google Scholar::
Simon SA, Collins TK, Kauffman GL, McNab WH, Ulrey CJ (2005). Ecological zones in the Southern Appalachians: first approximation. Research paper SRS-41, Southern Research Station, USDA Forest Service, Asheville, NC, USA, pp. 41.
::Online::Google Scholar::
Spaulding HL, Rieske LK (2010). The aftermath of an invasion: structure and composition of Central Appalachian hemlock forests following establishment of hemlock woolly adelgid, Adelges tsugae. Biological Invasions 12: 3135-3143.
::CrossRef::Google Scholar::
Stadler B, Müller T, Orwig D (2006). The ecology of energy and nutrient fluxes in hemlock forests invaded by hemlock woolly adelgid. Ecology 87 (7): 1792-1804.
::CrossRef::Google Scholar::
Thomlinson JR, Bolstad PV, Cohen WB (1999). Coordinating methodologies for scaling land cover classifications from site-specific to global: steps toward validating global map products. Remote Sensing of Environment 70 (1): 16-28.
::CrossRef::Google Scholar::
Tooke TR, Coops NC, Goodwin NR, Voogt JA (2009). Extracting urban vegetation characteristics using spectral mixture analysis and decision tree classifications. Remote Sensing of Environment 113 (2): 398-407.
::CrossRef::Google Scholar::
Trotter RT, Morin RS, Oswalt SN, Liebhold A (2013). Changes in the regional abundance of hemlock associated with the invasion of hemlock woolly adelgid (Adelges tsugae Annand). Biological invasions 15 (12): 2667-2679.
::CrossRef::Google Scholar::
USGS (2012). EarthExplorer. Web site.
::Online::Google Scholar::
Vapnik V (1995). The nature of statistical learning theory. Springer-Verlag, New York, NY, USA, pp. 189.
::Google Scholar::
Vastaranta M, Kantola T, Lyytikäinen-Saarenmaa P, Holopainen M, Kankare V, Wulder MA, Hyyppä J, Hyyppä H (2013). Area-based mapping of defoliation of Scots pine stands using airborne scanning LiDAR. Remote Sensing 5 (3): 1220-1234.
::CrossRef::Google Scholar::
Ward JS, Montgomery ME, Cheah CJ, Onken BP, Cowles RS (2004). Eastern hemlock forests: guidelines to minimize the impacts of hemlock woolly adelgid. Northeastern Area State and Private Forestry, USDA Forest Service, Morgantown, WV, USA, pp. 1-27.
::Google Scholar::
Webster JR, Morkeski K, Wojculewski CA, Niederlehner BR, Benfield EF, Elliott KJ (2012). Effects of hemlock mortality on streams in the southern Appalachian Mountains. The American Midland Naturalist 168 (1): 112-131.
::CrossRef::Google Scholar::
Wimberly MC, Reilly MJ (2007). Assessment of fire severity and species diversity in the Southern Appalachians using Landsat TM and ETM+ imagery. Remote Sensing of Environment 108 (2): 189-197.
::CrossRef::Google Scholar::
Wulder MA, Dymond CC, White JC, Leckie DG, Carroll AL (2006). Surveying mountain pine beetle damage of forests: a review of remote sensing opportunities. Forest Ecology and Management 221 (1): 27-41.
::CrossRef::Google Scholar::
Wulder MA, White JC, Coggins S, Ortlepp SM, Coops NC, Heath J, Mora B (2012). Digital high spatial resolution aerial imagery to support forest health monitoring: the mountain pine beetle context. Journal of Applied Remote Sensing 6 (1): 062527-1.
::CrossRef::Google Scholar::
Young RF, Shields KS, Berlyn GP (1995). Hemlock woolly adelgid (Homoptera: Adelgidae): stylet bundle insertion and feeding sites. Annals of the Entomological Society of America 88 (6): 827-835.
::CrossRef::Google Scholar::
Young JA, Morton DD (2002). Modeling landscape-level impacts of HWA in Shenandoah National Park. In: Proceedings of the “Hemlock Woolly Adelgid in the Eastern United States Symposium” (Onken B, Reardon R, Lashomb J eds). East Brunswick (NJ, USA) 5-7 Feb 2002. Agricultural Experiment Station, Rutgers University, New Brunswick, NJ, USA, pp. 73-85. -
::Online::Google Scholar::


Paper Contents

Paper Sections

Paper Figures

Paper Tables



Kantola T, Lyytikäinen-Saarenmaa P, Coulson RN, Holopainen M, Tchakerian MD, Streett DA (2016).
Development of monitoring methods for Hemlock Woolly Adelgid induced tree mortality within a Southern Appalachian landscape with inhibited access
iForest - Biogeosciences and Forestry 9: 178-186. - doi: 10.3832/ifor1712-008
First Previous Next Last
© iForest

Download Reference

Paper ID# ifor1712-008
Title Development of monitoring methods for Hemlock Woolly Adelgid induced tree mortality within a Southern Appalachian landscape with inhibited access
Authors Kantola T, Lyytikäinen-Saarenmaa P, Coulson RN, Holopainen M, Tchakerian MD, Streett DA
Close Download