Research Article, J Biodivers Manage Forestry Vol: 5 Issue: 1
Vegetation Cover Change in Glacier National Park Detected using 25 Years of Landsat Satellite Image Analysis
| Christopher Potter* |
| NASA Ames Research Center, Mail Stop 232-21, Moffett Field, CA, USA |
| Corresponding author : Christopher Potter NASA Ames Research Center, Mail Stop 232-21, Moffett Field, CA, USA Tel: 650-604-6164 E-mail: chris.potter@nasa.gov |
| Received: December 15, 2015 Accepted: January 22, 2016 Published: January 27, 2016 |
| Citation: Potter C (2016) Vegetation Cover Change in Glacier National Park Detected using 25 Years of Landsat Satellite Image Analysis. J Biodivers Manage Forestry 5:1. doi:10.4172/2327-4417.1000152 |
Abstract
Landsat satellite imagery at 30-m ground resolution was analyzed for Glacier National Park (GNP) Montana (United States of America) to determine when and where canopy green cover had changed significantly across the National Park area since the mid- 1980s. Image comparisons for the years 1984 and 2010 showed that consistent increases in the normalized difference vegetation index (NDVI) have been observed at several treeline areas, some of which have been monitored for increased tree cover for more than 75 years. Significant correlations between positive NDVI change and elevation were detected between 1500 and 1800 m elevation. Nonetheless, the greatest changes in NDVI over the past three decades in GNP have been caused by wildfires that burned more than 15% of the Park’s forested area over the past 10-15 years. As of 2010, less than 50% recovery of pre-fire canopy green cover was observed in most areas burned at high severity during 2003 and 2006. It is hypothesized that documented snowpack reductions and associated springtime temperature warming may be having detectable impacts on the slowing rate of post-wildfire regeneration rates of forests burned in GNP.
Keywords: Landsat; Glacier national park; Treeline; Wildfire; Climate change;Forest vegetation
Keywords |
|
| Landsat; Glacier national park; Treeline; Wildfire; Climate change; Forest vegetation | |
Introduction |
|
| The northern Rocky Mountains is a semi-arid region, supporting conifer forests at relatively high elevations, where cool temperatures and snow-dominated moisture sources can sustain tree growth during the short summer season [1]. Treeline vegetation change in Glacier National Park (GNP) Montana has been documented in numerous previous studies [2-4]. Roush et al. applied spatial analysis to digital photo-pairs at 12 sites within the alpine treeline ecotone of GNP to examine vegetation changes over more than 70 years [5]. Tree cover at the ecotones increased in 10 out of the 12 sites, with a mean increase of 60%, and new tree establishment occurred at all study sites. Increases in minimum summer temperatures and declines in early spring snow water equivalent (SWE) were linked too much of the increase in tree cover at these GNP ecotones. | |
| Late-20th century snowpack declines have been unprecedented in magnitude across the northern Rocky Mountains region when placed within the context of reconstructions over the past millennium [6]. Snowpack reductions and associated springtime temperature warming may be having important impacts on stream flow and water supplies across the western United States [7]. Projections of continued regional warming over the next 50-100 years have been associated with predictions for GNP and other protected wilderness areas of the northern Rocky Mountain region of longer fire seasons, larger annual area burned at high severity, higher basin-wide evapotranspiration (ET) fluxes from vegetation, and decreased river water flows [1,6]. | |
| The effectiveness of Landsat as a tool for monitoring landscapewide changes in the northern Rocky Mountain region has been demonstrated in several previous studies. For example, Jakubauskas and Price used Landsat images to estimate conifer forest age, tree height, and biomass. Landenburger et al. and Jewett et al. used Landsat images for whitebark pine mapping [8-10]. Annual vegetation ET flux in montane and sub-alpine vegetation communities of the western U. S. has been closely correlated with the satellite normalized difference vegetation index (NDVI) data from the NASA MODIS sensor [11]. | |
| The same NDVI from the NASA Landsat sensor has been shown to be a reliable index to monitor large-scale change in green vegetation cover and forest productivity, especially following disturbance in remote mountain areas of the western U. S. [12-20]. Results from these Landsat studies cited above have shown that canopy green leaf cover typically increases rapidly over the first five years following a stand-replacing disturbance, doubling in value by about 10 years after the disturbance, and then leveling off to approach pre-disturbance (mature) stand values by about 25-30 years after the disturbance event [21]. | |
| In the present study, Landsat satellite imagery at 30-m ground resolution was analyzed for the first time for all of GNP to determine the locations and timing of significant NDVI changes across the National Park area since the mid-1980s. The study design controlled for variations in yearly snow water accumulation, areas burned by wildfire, and elevation gradients. The study objectives that were intended to support National Park climate change assessments and planning were to: | |
| • Understand spatial variations in vegetation cover density changes over the past 25 years at treeline sites; | |
| • Compare vegetation cover changes in areas burned by wildfires over the past 25 years to nearby unburned areas of the National Park. | |
| It is worth noting that a continuous 25-year time-series analysis of Landsat NDVI image data carried out for national parks in California by Potter (2015a and 2015b) was mainly to reveal variations in vegetation cover density related to short-term (interannual) precipitation amounts and snow water accumulations [22,23]. Because the present study for GNP was designed to better understand major decades-long changes in vegetation cover density at treeline sites and wildfire areas, a paired-image comparison approach (photographic and satellite) was applied to overcome and control for continuous yearly variations in precipitation amounts and snow water accumulation. | |
Study Area Description |
|
| GNP is located in northwest Montana along the Canada-United States border and covers 4430 km2, ranging in elevations from 950 m along the Flathead River canyon to 3190 m on Mount Cleveland (Figure 1). Bordering the Park’s eastern boundary is the Blackfeet Indian Reservation, and the western and southern boundary is the Lewis and Clark National Forest and the Flathead National Forest. | |
 |
Figure 1: Map of the Glacier NP study area, with shaded elevation (left) and land cover types (right). The 1-km buffer zones around all of the treeline ecotone sites (Table 3) studied by Roush et al. (2007) [5] are outlined in red on the elevation map. Land cover color legend: Dark Green – Mature forest, Lighter Green – Successional forest, Yellow – Nonforest (shrubland and grassland), White – Ice and barren, Blue – Water, Red - Developed. |
| Before 1910, the GNP area was protected as a forest reserve [24]. GNP is still managed as a wilderness area, such that changes in vegetation cover would have been little affected by direct human activities, mainly logging or residential developments. Climate change, wildfire, and pest and insect outbreaks are the main remaining drivers of vegetation cover change in GNP, making it a logical choice for studies of natural moisture availability and NDVI trends in the northern Rocky Mountains region over the previous three decades. | |
| The Continental Divide of North America runs diagonally through GNP, creating a Pacific maritime climate to the west of the Divide, and a drier continental climate to the east. Winters in the area are long and cold, lasting from mid-November to late April. Summers are short and often dry, usually lasting from July through August. Bellaire et al. reported that air temperatures in GNP have increased between +0.5°C and +0.7°C over the past 30 years [25]. Roush et al. [5] reported that minimum summer temperatures at a station 70 km from GNP have increased about 1.5°C over the past century. | |
| Measurements of peak SWE levels at U. S. Department of Agriculture Snow Telemetry (SNOTEL) station at Flattop Mountain (Latitude 48.80oN, Longitude 113.86oW) within GNP confirm that the years 1992, 2001, and 2005 were among those with the greatest negative departures from mean snow pack levels recorded between 1980 to 2014 (Figure 2), while 1991, 1997, and 2011 showed the greatest positive departures from mean SWE levels [26]. Mean peak SWE was estimated at 116 cm, with a standard deviation of 26 cm. | |
|
Figure 2: Monthly SWE (1969 to 2014) recorded at the Flattop Mountain SNOTEL station within GNP. |
| The montane forest zone in the study area is situated between 1200 and 1800 m, and the subalpine forest zone is situated between 1800 and 2700 m, approaching treeline [27]. The vegetation mapping product from Hop et al. used photo interpretation and over 620 sampling plots within GNP to map 12 dominant National Vegetation Classification (NVC) System plant associations. GNP has more than 65% subalpine and montane forest cover, with alpine shrub, meadow, and talus cover at just over 16% of the total Park area (Figure 1 and Table 1). The major wildfire events that each burned more than 5000 ha within GNP over the past 25 years were summarized in Table 2 (from the National Monitoring Trends in Burn Severity, MTBS) [28]. | |
|
Table 1: Land cover class areas for GNP [24]. |
|
Table 2: Major wildfire events that each burned more than 5000 ha within GNP over the past 25 years. |
| Study sites selected by Roush et al. for repeat photographic analysis of GNP treeline ecotones (listed in Table 3) were located at the upper limit of tree cover and open subalpine meadows at slightly lower elevations (as shown in Figure 1). The dominant treeline species at all of these sites was subalpine fir (Abies lasiocarpa [Hook.] Nutt), with occasional Engelmann spruce (Picea engelmannii Parry ex Engelm.). | |
|
Table 3: Sites within the alpine treeline ecotone of GNP analyzed by Roush et al. [5] using ground-based repeat photography to detect changes in tree cover between (the base year) 1920 and 2003. |
Methods |
|
| Image processing | |
| Imagery from the Landsat Thematic Mapper (TM) 5 sensor was selected for the years 1984, 1989, 1994, 2003, and 2010 from the United States Geological Survey (USGS) Earth Explorer data portal (http://earthexplorer.usgs.gov/). These five years were similar in SWE, with peak monthly levels slightly below the 30-year average (from 1984 to 2014), as shown in Figure 2. Landsat image data from path/row 41/26 were consistently acquired for an anniversary window between July 15 and August 15 for all years, around the peak of the snow-free growing season in the northern Rocky Mountain region, to minimize variation caused by seasonal vegetation growth and sun angle differences [29]. | |
| All images used in this study were geometrically registered using terrain correction algorithms (Level 1T) applied by the USGS EROS Data Center. The Landsat Surface Reflectance Climate Data Record applied corrections to all the images used for top of atmosphere (TOA) reflectance, brightness temperature, and generated masks for clouds, cloud shadows, adjacent clouds, and surface water bodies [30]. | |
| Cloud-filtered NDVI (scaled from 0 to 1) was computed for all Landsat images as the differential reflectance between the red and near-infrared (NIR) portions of the spectrum by the equation: | |
| NDVI = (NIR–Red) / (NIR+Red) | |
| where NIR is the reflectance of wavelengths from 0.76 to 0.9 μm and Red is the reflectance from 0.63 to 0.69 μm. Low values of NDVI (near 0) indicate barren land cover whereas high values of NDVI (near 0.9) indicate dense canopy vegetation cover [8]. | |
| NDVI value differences were computed by raster image calculations for the dates selected. Paired image comparison approaches (whether photographic or satellite) are based the same assumptions and basic methodology, and carry similar potential for uniquely valuable results. The only major difference is that comparisons of two years of Landsat imagery can control for climate conditions in the selection dates of imagery, and have superior georeference capacity for the satellite image data values, compared to repeat photographic analysis. | |
| Spatial layers | |
| Elevation at 1 arc-second resolution was derived from the USGS National Elevation Dataset (NED). Vegetation cover types within GNP were determined by [24] based on aerial photographs. Areas within GNP that had burned in wildfires over the past century were delineated from the Fire History Polygons for the Northern Rockies 1889-2003, a multi-agency project designed to consistently map the burn severity classes (low, moderate, and high) and perimeters of fires across all lands of the United States [28,31]. | |
| Statistical analysis | |
| Tests of statistical significance between NDVI dates were carried out using the two-sample Kolmogorov-Smirnov (K-S) test, a nonparametric method that compares the cumulative distributions of two data sets [32]. The K-S test does not assume that data were sampled from Gaussian distributions (nor any other defined distributions), nor can its results be affected by changing data ranks or by numerical (e.g., logarithm) transformations. The K-S test reports the maximum difference between the two cumulative distributions, and calculates a p value from that difference and the group sample sizes. It tests the null hypothesis that both groups were sampled from populations with identical distributions according to different medians, variances, or outliers. If the K-S p value is small (i.e., p<0.05), it can be concluded that the two groups were sampled from populations with significantly different distributions. Within selected cover types of GNP, 1000 to 2000 point locations were randomly identified to test for differences between cumulative distributions of NDVI as a function of time, elevation, and vegetation cover types. | |
| Descriptive statistics of the distributions of sampled NDVI data sets were computed, including skewness and kurtosis [33]. Negative skewness indicates that the tail on the left side of the probability density function is longer than the tail on the right side, with values between -3 and +3 indicative of a normal distribution. Kurtosis is a descriptor of the shape of a distribution, with values greater than +3 indicative of a heavier tail than a normal distribution. | |
Results |
|
| NDVI changes at selected treeline sites | |
| A comparison of change over time in mean yearly NDVI from 1984 to 2010 for selected treeline sites described by Roush et al. [5] showed a gradual increase in NDVI between 1984 and 2010 at all of these sites (Figure 3). The period of most rapid increase in NDVI was generally observed between 1984 and 1994, followed by a slower rate of increase between 1994 and 2010. | |
|
Figure 3: Change over time in mean NDVI in GNP from 1984 to 2010 for selected treeline sites described by Roush et al. [5]; N >500 randomly selected points within each treeline 1-km buffer zone. Error bars for each year indicate 2 standard errors (2SE) of the mean NDVI per site. |
| The statistical comparison of NDVI differences between 2010 and 1984 sampled at the treeline sites listed by Roush et al. [5] showed that NDVI increased significantly (K-S p<0.01) at four of the nine sites where cloud-free Landsat data was available (Table 4). Site numbers 10-12 from Roush et al. [5] were completely obscured by clouds in either 1984 or 2010 images and hence could not be included in the NDVI comparison. Although positive NDVI differences between 2010 and 1984 at four of the nine cloud-free treeline sites were not statistically significant, they showed no overall tendency for declines in canopy green cover within the 1-km treeline buffer zone sampled. Only 3.7 ha around one (i.e., No. 6 - Granite Park) site, out of more than 11,300 ha sampled within these nine treeline buffer zones, had been affected by fires over the past century. | |
|
Table 4: Kolmogorov-Smirnov test results for the cumulative distributions of NDVI 1984 and NDVI 2010 for nine sites studied by Roush et al. [5] within the alpine treeline ecotone of GNP (N > 500 points sampled with cloud-free NDVI within each 1-km site buffer zone). |
| Maps of NDVI differences for the four treeline sties detected with significant positive canopy greening from 1984 and 2010 showed the extent of changes between +0.1 to +0.2 NDVI units per pixel (Figure 4). There were numerous clusters of pixel areas covering 0.5 km in length where NDVI was seen to increase by more than 0.1 unit since 1984. These clusters of pixels where the highest levels of positive NDVI change were almost all located along the margins of areas mapped by Hop et al. [24] as transitioning between sub-alpine forest and alpine scrub cover. | |
|
Figure 4: Changes in NDVI in GNP from 1984 to 2010 for selected treeline sites described by Roush et al, [5]. Shaded elevation is shown as the base map wherever cloud, ice, and water cover was detected and removed from the NDVI difference results. 1-km buffer zones were outlined in the circular polygons. |
| Only two (Nos. 2 Mokawantas Junction and 3 Haystack) of the nine Roush et al. [5] treeline sites sampled showed a notable positive NDVI change with elevation (Figure 5). In both cases, the significant correlation (R2>0.3; p< 0.05) between NDVI change and elevation was detected between 1500 and 1800 m. The strongest correlation between NDVI change and elevation was observed for the No. 3 Haystack site (Figure 4), where Roush et al. [5] measured a 6% increase in tree cover over the past 75 years. This treeline site had the lowest elevation range within the nine 1-km buffer areas sampled, with maximum elevation at around 1926 m (Table 3), whereas all the other treeline sites sampled had maximum elevation above 2200 m, and the majority of these sites had maximum elevation above 2500 m. | |
|
Figure 5: Correlation between NDVI change (2010 -1984) and elevation for two treeline sites described by Roush et al. [5]; N>660 randomly selected points within each treeline 1-km buffer zone. |
| NDVI changes across GNP | |
| The NDVI differences between 2010 and 1984 for all of GNP were dominated by areas burned by wildfires since the year 2000 (Figure 6). The majority of these large burned areas were detected with the 25- year NDVI difference of -0.1 to -0.2 units. Most notably, during the summer of 2003, nine wildfires burned nearly 700 km2 in and adjacent to GNP [34]. These 2003 fires affected roughly 13% of GNP’s total area [35]. The Robert Fire, which was the largest of the 2003 wildfires, alone burned 213 km2 on the west side of the Park and burned with high severity (according to MTBS burn classes), particularly on slopes along the northern shores of Lake McDonald [28,34]. The Red Eagle Fire of 2006 that extended across the eastern boundary of the Park also showed extensive negative NDVI differences of between -0.1 to -0.2 units at elevations between 1500–1900 m. | |
|
Figure 6: Changes in NDVI for the GNP area from 1984 to 2010. Cloud, ice, and water cover was masked as white. Boundaries of wildfires that have burned within GNP since the year 2000 were outlined in dark gray and the GNP boundary outlined in green. |
| High resolution (1-m) aerial imagery of selected areas that had burned by fires since the year 2000 revealed several noteworthy patterns of apparent vegetation response to disturbance (Figure 7). In the Wedge Canyon Fire example, areas at treeline (above 2150 m) that did not burn since 2000 commonly increased in NDVI by > +0.1 difference units, whereas many areas that did burn in 2003 had not recovered more than 50% of their pre-fire canopy green cover and appeared to be sparsely vegetated, compared to unburned forest areas nearby. The same pattern of less than 50% recovery of pre-fire canopy green cover was observed in the Robert Fire 2003 example shown at 1060 m elevation on northern-most shoreline of Lake McDonald. The example from the Red Eagle Fire of 2006 revealed a widespread burned area with less than 40% recovery of pre-fire canopy green cover by 2010 and extremely sparse regeneration of vegetation cover above 1700 m elevation. | |
|
Figure 7: True-color aerial imagery (at left) from the Department of Agriculture’s National Agricultural Imagery Program (Davies et al., 2010), paired with changes in NDVI from 1984 to 2010 (at right) for selected vegetation areas burned by fires since the year 2000. Boundaries of wildfires were outlined in both image types [43]. |
| The distribution of NDVI difference (2010 - 1984) values for all 30-m resolution pixels within the GNP boundary (Figure 8) showed that the areas burned by wildfires since 2000 had skewed the overall distribution slightly toward the negative extremes of <-0.1 NDVI difference units. Tests of normality for the distributions resulted in skewness of -2.5 and kurtosis of 6.5 for all 30-m pixels in GNP, skewness of -0.75 and kurtosis of 2.1 for randomly sampled NDVI difference values within areas burned since 2000, and skewness of -2.4 and kurtosis of 7.6 for areas that were unburned since 2000. Statistical comparison of NDVI differences (N=2000 random samples) between 2010 and 1984 showed that areas burned since the year 2000 were significantly more negative (K-S p < 0.001) than areas that were unburned since 2000, with a maximal difference of 0.42 NDVI units. | |
|
Figure 8: Histograms of change in NDVI for GNP from 1984 to 2010. |
| Points of interest within GNP | |
| A closer examination of the NDVI differences between 2010 and 1984 across all of GNP (Figure 6) showed several points of interest where NDVI had increased by at least +0.1 difference units over the 25 years of analysis. For example, areas that were burned by wildfires in both 1988 and 1994 in the northwestern section of GNP, all below 1600 m elevation, showed this level of strong recovery of green canopy cover by 2010. | |
| Several areas near treeline also showed widespread increases of > +0.1 NDVI difference units over the 25 years of analysis. Notable examples were in the headwaters of the Waterton River drainage above 1950 m elevation, along the nearby Livingstone Ridge to the south, and in the Yellow Mountain treeline area between 1970 and 2050 m elevation (all identified in Figure 6). | |
| The Belle River drainage in the northeast corner of GNP (Figure 6) showed among the highest consistent increases in (unburned) forest canopy cover across the entire park study area. NDVI differences between 2010 and 1984 were highest (>+0.2 units) along the river bottoms at around 1600 m elevation and remained high (>+0.1 units) up-slope to around 1780 m elevation. | |
| Numerous examples of a negative NDVI differences (<-0.1 units) between 2010 and 1984 were observed along the margins of sub-alpine lakes (such as Bowman and Logging Lakes and Lake Sherburne), located mainly in the northern sections of GNP. None of these examples had a record of wildfire around the water edges over the past century, leaving no apparent explanation for these changes in NDVI other than lake water level fluctuations. | |
Discussion |
|
| The Landsat satellite NDVI has been shown to be an effective tool to monitor large-scale change in forest productivity, especially following disturbance in remote mountain areas of the western U. S. [8,21-23]. Advantages of NDVI for the purpose of vegetation monitoring have been cited in its mathematical simplicity and ease of comparability across numerous multi-spectral remote sensing platforms [36]. Paired image comparisons of Landsat NDVI over 25- 30 years can add spatial details and valuable spectral information to historical photo-pair analysis that date back over seven decades [5]. | |
| To address the primary research objectives of better understanding spatial variations in vegetation cover density changes at treeline sites and in areas burned by wildfires over the past 25 in the National Park, the results of this study support the conclusion that the greatest changes in NDVI over the past three decades in GNP have been caused by wildfires that burned over the past 10-15 years within the Park. In contrast, significant increases in NDVI were detected at several treeline locations in GNP, as verified by the comparisons with long-term photo-interpretation results previously published by Roush et al. [5]. | |
| Specifically, the No. 3. Haystack treeline site was described by Roush et al. [5] as a valley area likely having had a greater potential for excessive snow drifting, which has been negatively correlated with new tree establishment by limiting the growing season length. Roush et al. [5] also made the point that snow accumulation can benefit established trees by protecting them from winter desiccation and injury. This treeline site was reported with one of the lowest rates of tree loss (4%) from the 75 years of photo-interpretation among all the sites studied by Roush et al, [5]. These ecological effects of relatively low-elevation valley topography resulted in the highest overall increase in NDVI from 1984 to 2010 among the treeline sites listed in Table 3. | |
| The No. 6 Granite Park treeline site was described by Roush et al. [5] as notable because 61% of the area changed over 75 years, largely as infilling (41%). This site also showed a relatively low level of new tree establishment (14%). Tree growth proceeded outward from the lower center of the site, where tree cover was historically the densest. The NDVI difference map result for site (Figure 4) is consistent with that spatial pattern of infilling and the 1-km buffer zone was among the highest overall for NDVI increase from 1984 to 2010 among the treeline sites listed in Table 3. | |
| For a more regional perspective on treeline dynamics, Romme and Turner [37] speculated that climate change could move the upper timberline in nearby Yellowstone National Park to 3300 m elevation or higher with projected warming trends over the next century. They further hypothesized that increased evapotranspiration, without compensating precipitation inputs or higher water use efficiency by plants, would likely increase drought stress in downstream river drainages. | |
| In sub-alpine zones of GNP where Landsat NDVI analysis did not detect an upward elevation shift in (unburned) canopy green cover near over the past 25 years, one hypothesis to explain this observation could be that whitebark pine (Pinus albicaulis), which retains snow and reduces erosion at high elevations while producing seeds that are an important food source for grizzly bears and other wildlife, is experiencing unprecedented mortality throughout its natural range [38]. A primary cause of this mortality in GNP is blister rust (Cronartium ribicola), an introduced pathogen that increases whitebark pine vulnerability to infestation by mountain pine beetle (Dendroctonus ponderosae) [39,40]. | |
| Declines in NDVI with dieback of whitebark pines are more likely to be detected during years of low moisture availably than they would be under more typical SWE levels of the past across northern Rocky mountain region [6,41]. Smith et al. [39] reported that treeline areas in GNP associated with high blister rust incidence were also found to be exposed to high wind speeds and were in close proximity to topographic depressions and wetlands. | |
| Climate change models have predicted average annual surface temperatures in the northern Rocky Mountain region to increase by 3.5°C by 2100, with the greatest levels of warming in the winter months [42]. This level of warming would likely promote widespread mountain pine beetle outbreaks that could eliminate whitebark pine from the region, with serious implications for grizzly bear food sources and the overall species composition of forests in protected national parks of the region [10,40]. | |
Conclusions |
|
| Landsat image comparisons over time for GNP showed that NDVI had consistently increased between the years 1984 and 2010 at several treeline sites. Significant correlations between positive NDVI change and elevation were commonly detected across the NP area between 1500 and 1800 m elevation. Nonetheless, the greatest changes in NDVI over the past three decades in the NP have been caused by wildfires that burned more than 15% of the Park’s forested area over the past 10-15 years. By 2010, NDVI analysis detected less than 50% recovery of pre-fire canopy green cover in most areas burned at high severity, in the years 2003 and 2006, suggesting that a relatively slow process of forest regeneration is underway in these disturbed areas [18]. | |
Acknowledgments |
|
| The author thanks Dan Fagre and Greg Pederson for helpful comments on an earlier version of this study. | |
References |
|
|
|
