Quaking Aspen in A High-Use Recreation Area: Challenges of People, Ungulates, and Sodium on Landscape Resilience

1. Introduction

Quaking aspen (Populus tremuloides Michx.; hereafter aspen) is the most widely distributed tree species in North America [1,2] and an icon of montane landscapes. Across the globe, aspen populations are in decline due to elevated climate warming, drought, fire suppression, land clearing, and browsing by ungulates [3,4,5,6,7,8]. In some cases, populations are in mild or rapid decline [9,10,11]. Specifically, in western North America, some stands are experiencing population declines [12,13], while others are shifting their range or increasing [14,15]. Indeed, both increased and decreased aspen coverage have been documented in the same landscape [14,16]. Hence, there is a need for site-specific and appropriate aspen management based on data-driven community assessments to address common threats, such as overbrowsing [17].

Aspen are of global environmental concern because they serve as a keystone species [18,19,20]. Stands of aspen support high biodiversity, representing an ecologically important habitat type for many plant and animal functional groups [20,21,22,23]. As a short-lived, clonal, broadleaf tree species with social value, aspen often serve as recreational destinations [20,24,25,26]. These stands are classified as stable or seral, meaning they will remain dominated by only aspen stems over time or they will progress toward canopies consisting of one-to-several conifer species, respectively [27,28]. This study will examine stable stands only, due to the make-up of the aspen community at our study site. Stable aspen stands are characterized by a complex, multi -height and -age stand structure due to continuous gap-phase patterns of mortality and regrowth [29,30,31,32]. Stable aspen stands with diverse height structures are more resilient to interspecific competition and mass insect attacks, stem or root pathogens, or large-scale blowdowns [28,30,33]. Achieving a diverse height structure of stable stand is inhibited by chronic ungulate herbivory [12,17].

Overbrowsing by ungulates is a well-documented threat faced by aspen forests in the western United States due to overpopulation, lack of predators, and lack of herd movement [12,17,34,35]. Overbrowsing prevents aspen regeneration from reaching the next stage class, reducing recruitment to inadequate levels for successful stand replacement [36,37]. Fencing in small areas of aspen is one tactic to reduce browsing [38,39,40,41,42,43], but it is not suggested as a long-term solution as fences restrict the movement of other large animals and it is aesthetically unfavorable [44].

Another potential obstacle for aspen stand replacement is soil contamination. The vehicle for such contamination is commonly water transported from adjacent locales via surface or subterranean vectors. Soil in forested areas may be impacted by mining, tainted ground-water, and other human development, with soil sodium (Na) levels being a main point of concern [45,46,47,48]. There is limited research on the long-term effects of Na and salinity on aspen growth and stand replacement, though it is known that other woody species are adversely affected by saline soils [49,50].

Recreation may also adversely affect forests, but relationships between recreation activities and aspen population dynamics has yet to be studied. Though studies examining aspen have been conducted in areas that are used for recreation, such as in Fish Lake National Forest and the Book Cliffs region in Utah, none directly have considered how recreation impacts the health of the forest or how the health of the forest impacts users [51,52]. Outdoor recreation often has direct and indirect impacts on ecosystem functions, processes, aesthetics, wildlife and, in turn, the ecological integrity of areas as recreational resources [53,54,55]. The ecological impacts of recreation have the potential to be both beneficial and harmful to the landscape, as in the case of wildlife disturbance [56]. Similarly, management actions taken to benefit a landscape can affect visitor experiences.

Here, we explored aspen status in a popular recreation area and the elements contributing to those conditions. Specifically, we investigated how overall condition, geographic variance, soil chemistry, and key indicators affected long-term aspen health alongside putative human development and herbivory impacts. The objectives of this study were to: 1) assess the current condition of a specific aspen-recreational ecosystem, 2) assess how forest characteristics vary across this system, 3) determine key environmental indicators and their interactions, and 4) propose appropriate site-specific management actions as viable options for managers to sustain aspen forests and recreation opportunities. This study serves as a reference for research and management in aspen communities, as well as other vegetation types, in high-use recreation areas that also may be impacted by adjacent development.

2. Methods

2.1. Study Site

This study took place in two adjacent high-use recreation areas in Summit County, Utah: Right Turn Sage (RTS) and Run-A-Muk (RAM), divided by the Olympic Parkway Road. RTS (453,016.04E 4,507,288.81N) is a trail area designated for hiking, mountain biking, and cross-country skiing with over 16 km of trails, while RAM (453,195.93E 4,507,142.45N) is an off-leash dog area (Figure 1). Unlike typical dog parks, RAM offers over five kilometers of well-maintained trail and 17.4 ha of enclosed forested and sage-covered hills. Both areas are managed by the Snyderville Basin Special Recreation District (hereafter Basin Recreation). According to the trail counters [57] maintained by Basin Recreation at each of the three entrances to RTS and RAM, in 2023 there were an average of 584 visitors per day at RAM and 69 at RTS, for a total of 623 daily visitors.

The RTS/RAM trail areas range from 1,963.8 m to 2,109.7 m in elevation, with an average of 2,036.8 m, and slopes ranging from 0.2% to 35.4% trending west/southwest. It covers 1,429.4 ha¹, with an estimated 11.6 ha¹ (~0.8%) of this area dominated by aspen forests. The surrounding, steeper hillsides are dominated by Gambel oak (Quercus gambelii Nutt.) and other tall shrubs. The area surrounding the aspen and oak forests consists of shrubs, grasses, and forbs. Other tree species are lightly scattered throughout the forest, including bigtooth maple (Acer grandidentatum Nutt.), Rocky Mountain juniper (Juniperus scopulorum Sarg.), subalpine fir (Abies lasiocarpa Nutt.), and white fir (Abies concolor Lindl.). Large herbivores in RTS/RAM include North American elk (Cervus elaphus L.), mule deer (Odocoileus hemionus Raf.), and moose (Alces alces L.). Though the study area was likely affected by historic livestock grazing, there are no recent uses for domestic animals in this area. From 1991 to 2020, the annual average precipitation for Snyderville, Utah is 44.34 cm mostly in the form of snow [58]. The average annual temperature was 6.78°C in neighboring Park City, Utah from 2006 to 2020 [58].

Utah Olympic Park, built for the 2002 Winter Olympics, is upslope and adjacent to the study area boundary and 1.6 km up the road from the RTS/RAM trail area parking lots. In August of 2019, residential housing for Utah Olympic Park was constructed, which altered the perennial and ephemeral path of the streams that flow through RTS/RAM, with the path now descending from a vehicle parking lot. Additional construction for the park is ongoing and may affect run-off content and patterns within the study area.

2.2. Field Methods

We collected data at RTS/RAM during July and August of 2022 using established aspen landscape survey methods [12]. To locate survey plots, we used ArcPro^® [59] to overlay a 30 m grid on the RTS/RAM boundary and randomly sub-selected 45 sample points (plots) that intersect a digital aspen cover layer. The aspen cover layer was created in ArcPro^® for this project using satellite imagery from the National Agriculture Imagery Program provided by the Utah Geospatial Resource Center [60]. At each plot, the environmental variables we recorded are sample area (RTS or RAM), number of vertical stand layers, a visual estimate of stand condition, understory cover, and percent aspen canopy cover.

All measures were assumed to represent conditions within approximately 0.5 ha¹ (radius ~40 m) surrounding the plot location center point . We identified vertical stand layers by looking horizontally through the forest and counting clearly distinguishable aspen layers across the entire site (i.e., short juveniles, tall juveniles, intermediate sub-canopy, and mature canopy stems). The presence of small numbers or patches of regeneration or recruitment did not a priori constitute an easily distinguishable “layer.” Where layers could not be determined due to continuous vertical stand structure, we recorded four layers - the maximum value.

Stand condition was determined based on a qualitative visual estimation of plot conditions developed as a time-saving method for forest assessment [61]. We used a standardized protocol for visually assessing tree damage, mortality, aspen layers, and overall browse impact to categorize aspen forest conditions as “Poor,” “Moderate,” or “Good” (Table 1). The Moderate category encompasses a much greater range of conditions, whereas Poor and Good groups are defined by proportional extremes.

Percent aspen canopy cover was measured along two 24 m transects using a densitometer. The transects ran east to west, were parallel, and 4 m apart. Aspen canopy cover was recorded at 1 m intervals along each transect, only counting live branches above 2 m height and live branches, resulting in fifty canopy cover readings from both transects that were then calculated into an overall stand percentage.

A fecal count of pellet groups was conducted for deer and elk along the same transects in 2 m belts to determine ungulate presence indirectly and approximate the level of aspen browse. No livestock are permitted within the study area, so these fecal types were not tallied. Pellet groups were defined as three or more pellets from the same defecation [62].

Stand structure is generally classified into three height stages: regeneration, recruitment, and mature overstory [12]. In this study, an individual aspen “tree” will be referred to as a “stem” - one stem of a single clone that is connected to the larger organism via an underground root network. Aspen stems within a 7 m radius of the center of each plot were tallied and classified as either regeneration (i.e., stems <2 m height), recruitment (i.e., stems >2 m height and <8 cm diameter breast height [DBH]), or mature (i.e., stems >2 m height, >8 cm DBH). Trees other than aspen were recorded noting their species, status (dead or alive), and DBH. Stems were considered separate individual stems if they forked below the soil litter. We categorized regeneration stems into three height classes (<0.5 m, 0.5-1 m, >1- 2 m), and leaders were examined for browsed buds and twigs to determine the percent browsed regeneration. For mature stems, we recorded the status (dead or alive), DBH, and damage. We recorded damage if >25% of bole circumference, crown, or root area was affected by stem cankers (fungal damage and/or resinosis), stem conks (>1 external conk), open wounds (exposed stem wood without decay from physical damage) root decay (not on the main bole; any visible fungi or exposed rotted roots), leaf damage (not on the main bole; spotted, rolled, eaten, or discoloration on >50% of leaves), or insect damage to bole (pin holes or obvious insect wounds or frass. We divided the number of mature stems with damage by the number of total mature trees to calculate the percent mature stem damage for each plot. All area-based data were summed, by plot, to the ha^-¹ level based on the area of sampling (i.e., either circular plot or rectangular belt transects).

Landscape position variables including elevation, slope, and aspect were derived from a 10 m digital elevation model (DEM) in ArcPro^® [59]. Aspect was transformed into a moisture index ranging from 0 (southwest aspect) to 1 (northeast aspect) to eliminate inaccuracies from averaging aspect values that straddle 0 degrees [63]. We also calculated the distance (m) from the road, hiking/biking trails, cross-country ski trails, and streams for each plot in ArcPro^®.

2.3. Soil Testing for Na, Ca, Mg, pH, SAR, and EC

A soil sample of ~300 g was collected 1 m north of each plot center from a depth of ~10-15 cm. Soil samples were air-dried, then a saturated paste extract (SPE) was prepared and used to measure the pH and electrical conductivity (EC). The SPE was analyzed by atomic absorption spectroscopy for calcium (Ca) and magnesium (Mg) and by flame emission spectroscopy for sodium (Na) using a 1:10 dilution [64]. The sodium adsorption ratio (SAR) was calculated to assess sodium hazard in soil solution [65]. SAR is the proportion of water-soluble Na to Ca plus Mg in the soil and is calculated with the following equation [66]:

$$SAR=\frac{Na}{\sqrt{\frac{Ca+Mg}{2}}}$$

(1)

The SPEs that yielded high results for Na (>20 mg/kg of soil) were further tested for pH, EC, Na, and 22 additional metals to ensure they were not contaminated during the original testing process [65,67]. Additionally, the original soil sample from the field were also further tested for pH, EC, Na, and 22 additional metals for SPE’s that yielded abnormally high Na levels (>100 mg/kg of soil) to ensure there was no contamination during the preparation of the SPE that would have causes high Na levels [65,67].

2.4. Data Analysis

We analyzed the data using R version 4.2.2 [68] and PC-ORD [69]. We considered results to be statistically significant at the 95% confidence interval (p <0.05; 69–71). Sampling values for all variables were compiled at the plot level for analysis. We divided plots located in the RTS portion into two elevation-based groups to provide sampling parity, resulting in three (overall) plot zones, as well as acknowledging an apparent natural division in water allocation (Figure 1). First, we ran preliminary descriptive statistics on plot-level data to identify outliers, test variance, and explore initial correlations in indicator pairings. We used non-parametric approaches due to the variation of data type and large numbers of zero values [70,72]. We mapped findings using ArcPro^® [59] for visual display and analysis. Analytical efforts for this work were of an exploratory nature [70]. We sought to determine the most important gradients among a suite of environmental variables (i.e., data reduction) that could serve as key metrics, or "indicators," of study area conditions.

Second, we used non-metric multidimensional scaling (NMDS) to explore correlations within this study's matrix of 45 plots and 26 response variables to highlight key variables [73]. The wide variation in data types/plot-level variables (e.g., stem counts, percent canopy cover, geometric location, elevation) requires a flexible and defensible analytical approach such as NMDS [74]. An initial outlier analysis was performed to check for data anomalies [74]. We initiated the ordination with a random start number upon 500 runs of the actual data set using the Sørensen distance measure. The final NMDS solution assessed dimensionality by plotting stress as a function of the number of dimensions or axes. Where two consecutive dimensions were <5 points of stress apart the lower dimension was selected as our optimum solution [70].

Third, multi-response permutation procedures (MRPP)—a non-parametric equivalent to MANOVA used for describing within-group agreement of variables in contrast to a priori data groups [75]—were used to test group difference for stand condition, sites, and geographic zones. We used R-package vegan to run the MRPP [76]. The Sørensen distance measure was used because it is less susceptible to bias due to outliers and zero values [74]. MRPP produces a T-value indicating the degree of difference between treatment group pairs, an A-value which is the chance-corrected within-group agreement (effect size), as well as a p-value establishing a level of test significance [70]. We tested two group configurations: sampling zones and Na levels categorized into low (<10 mg/soil kg), medium(10-20 mg/soil kg), and high(>20 mg/soil kg) Na content groups as recommended [77].

Lastly, we used two non-parametric tests to address indicators individually by significant groups found in MRPP testing: the Wilcoxon-Mann-Whitney U and Kruskal-Wallace tests [72]. We used the two-sided Wilcoxon-Mann-Whitney U test to evaluate field variables for differences between plots with adequate recruitment (>1200 stems ha^-1), sustainable browse levels (<20%) [78], and distance to the road or a stream (<5m and <10 m).

We employed the Kruskal-Wallace test, a non-parametric equivalent to analysis of variance, to test group differences between plots in different zones, with different stand layers, with different Na levels, and with different stand condition rankings. Additionally, the Kruskal-Wallace test served as the main method for evaluating the effectiveness and precision of the visual stand condition ranking (see Table 1 for ranking criteria). Stand condition is an amalgamation of recruitment ha^-1, dead basal area, browse, percent mature stem damage, and stand layers, so we used the Kruskal-Wallace test if these variables varied significantly between plots ranked as having “Poor,” “'Moderate,” and “Good” stand conditions (Table 1). Confirmation or rejection of qualitative stand ratings, based on actual field measures, is intended to elucidate usefulness of a quick visual rating as a reasonable substitute for more time-intensive field measures.

3. Results

3.1. Descriptive Statistics of Site Conditions

We examined the 45 plots across RTS/RAM by three zones to further understand geographic discrepancies, if any, within the study area: Zone 1, Zone 2, and Zone 3 (Figure 1). Table 2 displays mean and standard deviations for 29 variables by zone as well as for the overall study area. Notably, regeneration ha^-1 had a mean four to eight times greater than mean recruitment ha^-1 across the study site. Also, the mean recruitment ha^-1 does not meet the minimum stand placement of 1,200 ha^-1 [61,79], but the ratio of recruitment ha^-1 compared to mature ha^-1 indicated that there is enough recruitment to replace the existing overstory. Lastly, there are high amounts of Na [77] and percent browse [61].

Results from the soil testing yielded unusually high results for Na at plot 15, and high results for six additional plots. Figure 2 displays Na levels at each plot, with red and black plots having the highest levels of Na and green plots having the lowest. Visually, sample locations with high values appear closer to streams. This assertion was tested further below (see 3.4; Figure 6).

3.2. Ordination Seeking Key Site Indicators

The NMDS ordination produced a three-dimensional solution on a matrix of 45 plots by 26 variables, including Zones, for the 2022 dataset with a final stress value of 11.65 with an instability of 0.009. Stability was reached at 39 iterations from a maximum of 400 runs of our "real" dataset. A Monte Carlo test of 400 random data runs versus the real data set verified a significant NMDS outcome (p = 0.0349). No plots were eliminated in the outlier analysis, but the variable regeneration was eliminated in the outlier analysis (>2 SD, outlier analysis). Elevation was also removed due to insignificant variation in elevation across the site (Table 2). The recruitment-mature ha^-1 ratio and live-dead basal area ratio variables were also eliminated from the ordination because it contained null values.

The three-axis result described 88.1% of ordination variance (Axis 1: r² = 0.405; Axis 2: r² = 0.319; Axis 3: r² = 0.158, orthogonality Axis 1 and 2 = 96.5%). Although the significance of each axis does not meet our 95% confidence interval standard (Axis 1 = 0.09, Axis 2 = 0.04, Axis 3 = 0.03), we concluded that the significance is adequate, but should factor into our conclusions. Table 3 presents all environmental variables with their Pearson’s coefficient (r) values by axes identified in NMDS. The NMDS is displayed as one joint plot with Axis 1 and Axis 2 where vectors with > +/- 0.5 Pearson’s coefficient (r) value are displayed as an overlay in Figure 3. The length and direction of vectors in Figure 3 correspond to variable strength and relationship to the two-dimensional plot-data space. Axis 1 describes a gradient of recruitment abundance negatively correlated to deer scat presence. Axis 2 displays a gradient of canopy cover, mature ha^-¹, total basal area, and live basal area. On Axis 3, three variables have high Pearson’s coefficient (r) values, elk scat ha^-¹ (r = 0.518), deer scat ha^-¹ (r = 0.685), and Na (r = 0.32). Several variables do not meet our threshold for Figure 3 joint plot display, yet they provide support in their statistical relationships to other variables along dominant ecological gradients. In the joint plots, Zone 1 has less recruitment and poor overstory growth while Zone 2 trends in opposition to Zone 1. Zone 3 has no dominant trend.

3.3. Between-Group Difference for Zones and Sodium

We examined several categorical variables for group differences in plots using multi-response permutation procedures (MRPP). The most telling group differences were found by zone and soil Na. When comparing zones, results show that there is greater within group agreement (validation) than between groups for all three pairwise alignments (Table 4). Although comparisons are significant (<0.05) for Zone 1 vs 2 (p = 0.025) and Zone 2 vs 3 (p = 0.008), low T-values (Zone 1 vs Zone 2 = -2.49; Zone 2 vs Zone 3 = -3.4) indicate modest between group differences. Similarly, when comparing plots with low (<10 mg/soil kg), medium (10-20 mg/soil kg), and high (>20 mg/soil kg) Na categories, the MRPP results show that there is stronger within group agreement (validation) than between groups with low and high Na levels (Table 4). Comparisons are significant only between low and high Na categories (p = 0.02), but the low T-values (low vs medium = -0.83; low vs high = -2.7; medium vs high = -0.45) indicate weak between group differences. In other words, there are significant differences between the three zones and low vs high soil sodium levels, but those differences are not large.

3.4. What Indicators Cause between-Group Differences?

Figures 2.4, 2.5, and 2.6 present an array of non-parametric test results—Kruskal-Wallis test for three-way and Mann–Whitney U test for two-way comparisons—to discern the status of recruitment ha^-¹, percent canopy cover, and stand layers in each Zone (Figure 2.4), ungulate presence through percent browse, elk scat ha^-¹, and dear scat ha^-¹ in each Zone (Figure 5), and Na and SAR levels in relation to distance to stream (Figure 6).

Zone 1, the area of initial concern for managers due to apparent mortality, has the lowest mean for canopy cover, stand layers, recruitment, mature, live basal area, and dead basal area (Table 2). Zone 2 has the highest mean for elk scat, deer scat, and browse level (Table 2). Zone 3 has the lowest mean for elk scat, deer scat, and browse, and the highest mean for canopy cover, stand layers, recruitment, mature, and dead basal area (Table 2). Of these, statistically significant zonal differences are only found when comparing stand layers (X ² = 7.01, p = 0.036; Figure 4b), live basal area (X ² = 13.732, p = 0.001; Figure 4c), and deer scat (X ² = 11.92, p = 0.003; Figure 5c). There are trends that are not significant but interesting include the higher elk presence in Zone 2, low recruitment in Zone 1, and low ungulate visitation in Zone 3.

Figure 4. Box plots of nonparametric tests describe study site patterns of regeneration, recruitment, and browsing impacts by Zones. Analyses graphics are as follows: (a) difference in recruitment ha^-¹ for all study plots; (b) significant difference in stand layers for all study plots; (c) significant difference in live basal area for all study plots. Kruskal-Wallis test results are shown for differences between the three groups including chi-squared (X²) and significance values (p). The x-axis shows the Zone and the y-axis reports Wilcoxon mean scores. Output from Kruskal–Wallis test whiskers show minimum and maximum values, open circles are outliers, boxes represent 25%–75% data ranges, horizontal lines within boxes are medians, and open triangle symbols are mean.

Figure 4. Box plots of nonparametric tests describe study site patterns of regeneration, recruitment, and browsing impacts by Zones. Analyses graphics are as follows: (a) difference in recruitment ha^-¹ for all study plots; (b) significant difference in stand layers for all study plots; (c) significant difference in live basal area for all study plots. Kruskal-Wallis test results are shown for differences between the three groups including chi-squared (X²) and significance values (p). The x-axis shows the Zone and the y-axis reports Wilcoxon mean scores. Output from Kruskal–Wallis test whiskers show minimum and maximum values, open circles are outliers, boxes represent 25%–75% data ranges, horizontal lines within boxes are medians, and open triangle symbols are mean.

Figure 5. Box plots of nonparametric tests describe study site patterns of regeneration, recruitment, and browsing impacts by Zones. Analyses graphics are as follows: (a) difference in percent browsed for study plots; (b) difference in elk scat ha^-¹ for all study plots; (c) significant difference in deer scat ha^-¹ for all study plots. Kruskal-Wallis test results are shown for differences between the three groups including chi-squared (X²) and significance values (p). The x-axis shows the Zone and the y-axis reports Wilcoxon mean scores. Output from Kruskal–Wallis test whiskers show minimum and maximum values, open circles are outliers, boxes represent 25%–75% data ranges, horizontal lines within boxes are medians, and open triangle symbols are mean.

Figure 5. Box plots of nonparametric tests describe study site patterns of regeneration, recruitment, and browsing impacts by Zones. Analyses graphics are as follows: (a) difference in percent browsed for study plots; (b) difference in elk scat ha^-¹ for all study plots; (c) significant difference in deer scat ha^-¹ for all study plots. Kruskal-Wallis test results are shown for differences between the three groups including chi-squared (X²) and significance values (p). The x-axis shows the Zone and the y-axis reports Wilcoxon mean scores. Output from Kruskal–Wallis test whiskers show minimum and maximum values, open circles are outliers, boxes represent 25%–75% data ranges, horizontal lines within boxes are medians, and open triangle symbols are mean.

Figure 6. Box plots of nonparametric tests describe sodium (Na), sodium adsorption ratio (SAR), and meters from stream. Analyses graphics are as follows: (a) significant difference in meters from streams for study plots with Na levels of <10mg/soil kg, 10-20mg/soil kg, and >20mg/soil kg; (b) significant difference in levels of Na for study plots <10 m, 10-20 m, and >20 m from streams; and (c) significant difference in SAR for study plots <10 m, 10-20 m, and >20 m from streams. Kruskal-Wallis test results are shown for differences between the three groups including chi-squared (X²) and significance values (p). The x-axis shows the Zone and the y-axis reports Wilcoxon mean scores. Output from Kruskal–Wallis test whiskers show minimum and maximum values, open circles are outliers, boxes represent 25%–75% data ranges, horizontal lines within boxes are medians, and open triangle symbols are mean.

Figure 6. Box plots of nonparametric tests describe sodium (Na), sodium adsorption ratio (SAR), and meters from stream. Analyses graphics are as follows: (a) significant difference in meters from streams for study plots with Na levels of <10mg/soil kg, 10-20mg/soil kg, and >20mg/soil kg; (b) significant difference in levels of Na for study plots <10 m, 10-20 m, and >20 m from streams; and (c) significant difference in SAR for study plots <10 m, 10-20 m, and >20 m from streams. Kruskal-Wallis test results are shown for differences between the three groups including chi-squared (X²) and significance values (p). The x-axis shows the Zone and the y-axis reports Wilcoxon mean scores. Output from Kruskal–Wallis test whiskers show minimum and maximum values, open circles are outliers, boxes represent 25%–75% data ranges, horizontal lines within boxes are medians, and open triangle symbols are mean.

Levels of Na compared to distance from stream was significant in three different comparisons. First, between Na categories and distance from stream (X ² = 6.482, p = 0.014), with the Na categories being <10 Na mg/soil kg, 10-12 Na mg/soil kg, and >20 Na mg/soil kg (Figure 6a). Plots with >20 Na mg/soil kg are within 30 m of a stream, while 75% plots in the other two categories were between 30 m and 80 m from a stream. Second, between distance from stream categories and Na (X ² = 8.941, p = 0.011), with the distance categories being <10 m, 10-20 m, and >20 m (Figure 6b). Plots <10 m from a stream range between 2.38 and 132.65 Na mg/soil kg. Plots 10-20 m from a stream range between 2.31 and 47.9 Na mg/soil kg, with an outlier of 83.79 Na mg/soil kg. Plots >20 m from a stream were under 41.79 Na mg/soil kg but had a mean and median bellow 10. Lastly, between the same distance from stream categories and SAR (X ² = 6.626, p = 0.036; Figure 6c). Similar to plots with high Na, plots with an SAR between 0.20 and 8.15 are <10 m from a stream, while plots 10-20 m from a stream have an SAR of 0.19-2.37, and plots >20 m from a stream have an SAR of 0.17-4.28. In summary, there is a pattern of plots with high SAR and Na near streams. Plots with a high SAR have a high amount of Na compared to Mg and Ca, which could lead to the inability of water to penetrate soil.

The Kruskal-Wallace test was used to test the accuracy of stand condition rankings that were visually assessed in the field (Table 1). Test results show only one variable is statistically significant when comparing differences in plots with “Poor,” “Moderate,” and “Good” stand condition: recruitment ha^-1 (X ² = 6.3511, p = 0.042). The other variables that contribute to stand condition are not statistically significant: dead basal area (X ² = 3.591, p = 0.166), browse (X ² = 0.706, p = 0.702), percent mature stem damage (X ² = 1.664, p = 0.435), and the number of stand layers (X ² = 5.157, p = 0.076).

4. Discussion

Findings from the analysis indicate that ungulate presence (Figure 3) and soil contamination (Figure 1 and Figure 6) are negatively affecting aspen conditions at some localities, though much of the study area is relatively healthy. Browse and Na levels are variables that could explain the disconnect between regeneration stems and recruitment stems—i.e., there is substantial regeneration that is not making it to the next height class. Additionally, maintenance done for recreation may also be impacting these forests. Ungulate presence, Na levels, and recreation use vary by zone and impact aspen in unique, and potentially additive, ways.

4.1. How Does Aspen Condition Vary across the Study Area?

Recruitment stems are a result of the successful growth of regenerating stems into subcanopy trees and is often used to measure structural diversity and “escape” from the reach of browsers [12]. Low or absent recruitment strongly suggests a temporal pattern of decline by an inability to replace overstory mortality [12,30,35,36]. The most prominent result of ordination analysis here describes a clear inverse relationship between ungulate presence and successful aspen recruitment throughout our study area (Table 3; Figure 3). For the entire study site, mean recruitment does not exceed 1,200 stems ha^-¹—the amount of recruitment ha^-¹ considered the minimum for stand replacement (Table 2) [61,79]. However, the mean ratio of recruitment ha^-¹ versus mature ha^-¹ is greater than one for all three zones (Table 2), suggesting that there is at least a minimum level of juvenile stems for replacement should the overstory die-off [61]. This cursory assessment, though, is conservative and does not account for rapid mortality or attrition among these saplings prior to their reaching maturity or the geographic disparities we found within subsections of our sampling zones (i.e., individual plots with very poor recruitment). The amount of regeneration could be a result of gap-phase regeneration occurring because there are canopy openings allowing more sunlight, the mortality of mature aspen stimulating new growth, or the continuous stunting of regeneration stems due to overbrowsing [33,61,80]. The latter could explain why the level of recruitment ha^-¹ is not enough for adequate stand replacement despite apparently sufficient regeneration (Table 2).

Zones across the site were significantly different based on an amalgam of all measures taken (Table 3). Significant variables contributing to these differences include stand layers (Figure 4b), live basal area (Figure 4c), and deer scat (Figure 5c) being significant indicators, with Zone 1 having the least amount of stand layers and live basal area, and Zone 2 having the least amount of deer scat. Moreover, at least for Zone 1, most sample plots displayed negative patterns for successful recruitment and overall aspen growth (Figure 3). These results suggest that a further understanding of geographic variability—in recreational influences, herbivory patterns, and soil toxicity—within this study area is required to fully understand broad-scale impacts.

4.2. What Effect Does High-Use Recreation Have on Ungulate Presence and Aspen Conditions?

There is a clear pattern of less ungulate presence where dog walking is high at our site. In Zone 3, the off-leash dog park, the lowest mean elk scat, deer scat, and percent browse (Table 5). The smaller degree of ungulate deer scat, elk scat, and browse (Figure 5) could be due to higher recreation use and off-leash dogs or the scents of the dog urine and feces. Free-roaming dogs harass, compete with, and prey upon wildlife which alters the activity patterns of wildlife [53,54]. Visscher et al. [81] found that human impacts on deer use and behavior is greater than that of natural predators (i.e., canids, felids, and ursids), though the same was not true for elk. Other studies found that ungulates, specifically elk, have been documented to avoid areas where humans and dogs are present, though they return during unoccupied hours [82,83], which could be the case at our study site. In the current study, this means that recreators and their dogs could be deterring ungulates from overbrowsing on young aspen, having an indirect positive impact on aspen condition in a portion of this area.

Hunted populations of ungulates have more of a flight response to humans and dogs than unhunted ungulates [83]. Such a refugia affect for stalked wildlife in areas off limits to hunting can result in elevated impacts to browsed vegetation [51].This is especially true when ungulate and human populations in the area are high, as they are at our site, and habituation in heavy traffic areas may decrease avoidance behavior of ungulates [84]. Hunting is not permitted at our study site due to safety concerns of a recreation area, so the unhunted ungulates may be more habituated to human and dog presence in the area. In contrast, a recent study by Beirne et al. [85] shows that large herbivores in high human populated areas became more nocturnal as human activity increases. So it is possible that the high presence of humans and dogs in Zone 3 (Figure 5) is deterring deer for long enough to meaningfully reduce browsing. In any case, the relationship of recreationists with abundant off-leash dogs, their impact on ungulate movement and behavior, and in turn the change in the amount of browsing on aspen is a topic worth investigating further.

4.3. How Is Soil Contamination Impacting Aspen Communities?

At our study site, indicators of soil contamination occurs where aspen recruitment and mature stems are lowest (Figures 2.4 and 2.6). Zone 1 has the highest mean of SAR (1.07) and lowest mean of canopy cover, recruitment ha^-1, mature ha^-1, and live basal area (Table 2). SAR is the concentration of Na in soil measured as a ratio compared to magnesium (Mg) and calcium (Ca) [66,86]. A high ratio of Na compared to Mg and Ca, or a high SAR, is a concern for soil structure because it causes soil particles to disperse instead of clumping or aggregating [66,86]. Aggregated soils have larger pores, allowing water and root systems to move through the soil, whereas dispersion creates small pores that restrict the movement of water and root systems [65,86]. Dispersion restricts the movement of roots and water [86].

This is particularly apparent at plot 15 in Zone 1, a location along a stream with pooled water, creating a marshy area with broadleaf cattails (Typha latifolia L.). Plot 15 has the highest SAR (6.31), the highest Na (132.65 mg), and a low electrical conductivity (EC) (0.2 dS/m; Figure 2). The combination of an SAR >5 combined with an EC <2 dS/m results in ponding – regularly concentrating water on top of the soil [66,86]. We found a clear pattern of high Na near the streams at our study site. Figure 2 shows us that the seven highest Na levels (>20 mg) are near the waterways, indicated by warm colored dots and a black dot for plot 15. Additionally, Figure 6 show us that high Na and SAR levels occur most frequently with closer proximity to surface water.

The proximity to streams leads us to speculate about the upstream inputs into the water. The stream that flows through Zone 1 comes directly out of the parking lot for a housing development for the Utah Olympic Park, which was completed in 2019, before which the footprint of that land was undeveloped (see Figure 7). Additional Na could be added to the system from the construction of the housing unit, continued construction from the expansion of the Utah Olympic Park, or winter road salting intended to increase traction under icy conditions. Unfortunately, the source of additional Na in the system was beyond the scope of this study, although we are compelled to report that it appears to be negatively affecting aspen forests here. These results suggest further testing of soil contents and investigation into the origins of likely off-sight contamination is recommended.

EC in higher concentrations can reduce plant growth through toxicity effects. High levels of EC induce osmotic stress and limits the ability of roots to intake both water and nutrients, causing plants to wilt even when water is present [66,87,88,89]. Additionally, high levels of SAR and EC can limit root density, root length, and the ability of roots to intake both water and nutrients [87,88].

However, Lillies et al. [90] examined the growth of aspen in naturally saline sites, finding that aspen can grow in soil conditions previously considered unsuitable for forest vegetation (with an EC of 10 dS/m, and SAR of 13) with little evidence of nutritional toxicities or deficiencies. Long term effects of increased salinity included a significant decrease in basal area and overall growth over time. This matches our findings of Zone 1, the area of initial concern, having the lowest mean live basal area and mean dead basal area and the highest mean SAR compared to the other two zones (Table 2). Khasa et al. [91] looked at the effects of salinity on aspen seedlings in short-term greenhouse and field trials, but the results do not translate to aspen landscapes or show long-term effects. Vaario et al. [92] found that aspen experiencing stress from saline soils had lower leaf count on seedlings. Additionally, leaf necrosis significantly increased with increased salinity. However, this was a short-term study that did not examine long-term effects. Though the relationship specifically between aspen condition, Na, SAR, and EC has not been explored on a long-term basis, impenetrable soil caused by high SAR and low EC is likely important because it restricts the movement of water and roots, which is evident at our study site.

4.3. How Might Area Maintenance Impact the Aspen Landscape?

Recreation is shaping the ecology of our landscape in several ways. Above we posited that the high presence of humans and dogs in Zone 3 could deter ungulate browsing. We also suggest that construction or the application of road salt occurring upstream for recreation at Utah Olympic Park may be the source of sodium at our site. Maintenance done onsite for cross-country skiing is another factor worth considering. Managers mow paths that are ~6 m wide in Zones 1 and 2 roughly twice a year. Aspen regeneration is not exempt from this mowing, which prevents the establishment of any aspen in the designated cross-country ski paths. Repeat mowing may act to deplete carbohydrate reserves in the clonal root network, much like chronic herbivory [93]. Such denudation of energy reserves ultimately affects the ability of new suckering to grow and thrive. Twelve out of 34 plots in Zones 1 and 2 combined are within 1 m of or intersect the cross-country ski path (Table 1), which means that regeneration stimulated by dying trees within the plots are removed by mowing for at least a third of all plots in Zones 1 and 3. We are unsure of the degree of impact mowing has on the site as a whole as it was not a focus of this study and recommend further investigation.

5. Conclusion and Management Implications

We examined the condition of a unique, high-use aspen recreational landscape using established aspen landscape survey methods [12]. Across the study site, there is a lack of the vertical diversity expected of stable aspen stands for sustainable stand maintenance [29,30,31]. We suggest that management actions be paired with continued monitoring of the study site to understand system responses alongside ongoing soil testing with varying proximity to the streams. Temporary fencing of small areas is recommended to reduce browsing in Zone 1 and Zone 2 and to further understand the degree to which browsing is contributing to aspen conditions. Similarly, implementing a retention pond to filter out upstream pollutants is recommended to prevent further increases of Na into the system and to examine how soil SAR levels and EC change with the treatment. Managers could also consider permitting or encouraging users with off-leash dogs to visit RTS, where the apparent decline is occurring, for certain periods of time.

This aspen system is unique as it serves as a recreation destination for the local population, specifically for the off-leash dog area, and because of its location directly downstream from continued development. It should be clear from this and other works [29,51,94] that aspen vary considerably by topographic and geographic locale; thus, they should not be managed as one type. This systematic survey of a highly recreated landscape can serve as a reference when examining high-use and near-development aspen communities, as well as forested systems at-large.