Background

Meadow ecosystems have important ecological functions and support socioeconomic services, yet are subject to multiple stressors that can lead to rapid degradation. In the Sierra Nevada of the western USA, recreational pack stock (horses and mules) use in seasonally wet mountain meadows may lead to soil trampling and meadow degradation, especially when soil water content is high and vegetation is developing.

Methods

In order to improve the ability to predict meadow vulnerability to soil disturbance from pack stock use, we measured soil resistance (SR), which is an index of vulnerability to trampling disturbance, at two spatial scales using a stratified-random sampling design. We then compared SR to several soil and vegetation explanatory variables that were also measured at the two spatial scales: plant community type (local scale) and topographic gradient class (meadow scale).

Results

We found that local-scale differences in drivers of SR were contingent on the meadow scale, which is important because multiple spatial scale evaluation of ecological metrics provides a broader understanding of the potential controls on ecological processes than assessments conducted at a single spatial scale. We also found two contrasting explanatory models for drivers of SR at the local scale: (1) soil gravimetric water content effects on soil disaggregation and (2) soil bulk density and root mass influence on soil cohesion. Soil resistance was insufficient to sustain pack stock use without incurring soil deformation in wet plant communities, even when plant cover was maximal during a major drought.

Conclusions

Our study provides new information on seasonally wet meadow vulnerability to trampling by pack stock animals using multi-scale drivers of SR, including the contrasting roles of soil disaggregation, friction, and cohesion. Our work aims to inform meadow management efforts in the Sierra Nevada and herbaceous ecosystems in similar regions that are subject to seasonal soil saturation and livestock use.

Field sampling

We used a stratified random sampling design within meadow PCTs (representative of differing hydrologic regimes) to examine local-scale differences in SR and potential explanatory factors of SR along a soil moisture gradient from wet to dry within meadows. During the peak growing season in July 2013 (when plant cover is maximal), we sampled within 1-m² plots placed in four meadow PCTs (i.e., strata). We randomly located plot locations prior to fieldwork using ArcGIS 10.1 software, spacing plots more than 4 m apart from each other to ensure spatial independence of samples (Weixelman and Riegel 2012). Of the four PCTs selected for study, all four plant community types do not typically occur together in any one meadow. Thus, for each of the five study meadows chosen, we sampled up to three PCTs, taking care to ensure equal sample sizes for each.

We used meadow gradient class to examine meadow-scale differences in SR among meadows of differing topographic gradients, based on the hydrogeomorphic type (HGM) classification developed by Weixelman et al. (2011). Topographic gradient is important because soil moisture conditions (i.e., soil water content and surface water velocity) are in part a function of gradient. We examined three low gradient and two middle gradient meadows, with at least three PCTs sampled within each meadow (Fig. 2). To capture soil conditions within meadows, we sampled six randomly placed replicate plots per PCT during each of three sampling periods during peak growing season in July of 2013, resulting in 18 replicates total per PCT per week (total plots = 252, Fig. 2). We sampled soils in plots over 3 weeks due to the remoteness of study locations, where only a small number of samples could be carried out of the field weekly.

We oriented 1 m² quadrat plots north to south to avoid sampling bias and placed plots within each PCT representative of a soil moisture gradient from wet to dry within meadows. We required a minimum absolute vegetation cover of 25% for Mesic and the Xeric PCTs (CM, SK, DC), and a minimum cover of 15% for the hydric PCT (CV), which typically has less dense vegetation cover. Within the 1-m² sample plots, we used a nested sample-plot design of 15-cm and 5-cm radius circular plots (~ 700 cm² and ~ 80 cm², respectively) to measure SR, soil volumetric water content (VWC), and total vegetation cover (TVC). In the 1-m² and ~ 700-cm² plots, we ocularly estimated TVC, in order to more closely spatially pair SR data with VWC and TVC data measured. We ensured that the same observer performed ocular estimations of TVC throughout the study to avoid observer bias (Elzinga et al. 1998). In the ~ 80-cm² plots, we measured soil VWC and SR and collected undisturbed soil cores.

We measured SR at the center of nested circular plots, located 30 cm from the SW corner of each quadrat. To measure SR, we used a dynamic cone penetrometer (Synergy Resource Solutions Inc., Belgrade, MT, USA) based on the design from Herrick and Jones (2002). We measured soil VWC (0–12 cm depth) within a 5-cm radius to the northeast of penetrometer readings using a HydroSense II Time Domain Reflectometer (Campbell Scientific, Logan, UT, USA).

To evaluate potential drivers of SR, we took soil samples to assess soil physical properties that may covary with SR. We collected intact, undisturbed mineral soil core samples (5 cm diameter, 0–15 cm depth) using an AMS soil core sampler (AMS Inc., American Falls, ID, USA) within a 5-cm radius to the southwest of penetrometer SR readings. We extracted a total of 252 soil cores (one core per plot; Fig. 2). We enclosed intact soil cores in sealed 4-mm-thick polyethylene bags, kept cool in coolers with ice, and transported them to the laboratory at the University of California, Merced. At the laboratory, we stored soils at 4 °C for 6 months prior to processing, which occurred over a 1-month period.

Laboratory analyses

In the laboratory, soil was removed from the polycarbonate sleeves and sieved field-moist through a 4-mm mesh sieve, followed by sieving through a 2-mm sieve. Plant roots collected on the sieves were hand separated (mostly fine, i.e., < 2 mm diameter) and then oven-dried at 70 °C. We chose to sieve rather than elutriate samples for obtaining roots because sieving was the faster method. Coarse fragments (> 4 mm and 2–4 mm) were oven dried at 105 °C and weighed. We separated coarse fragments into two size classes to test whether these classes corresponded with differences in VWC and SR. Root content and coarse fragment contents (CF > 4 mm, CF 2–4 mm) were expressed as dry mass of material per dry mass of all material within the core multiplied by 100% (i.e., roots %). Root mass areal density (RMAD; dry root mass per unit ground area) was calculated by dividing the oven-dry root mass per core by the cross-sectional area of the core (19.6 cm²).

Soil gravimetric water content (GWC) was determined by drying subsamples (~ 5 g) of field-moist, sieved soil (< 2 mm) in an oven (105 °C for 48 h). We measured GWC in addition to field-measured volumetric water content (VWC) to determine if there were differences between the two. Soil bulk density (BD) was calculated as the oven-dry mass of soil (< 2 mm) per core volume. We also measured soil water holding capacity (WHC), modified from the procedure reported by Haubensak et al. (2002). Briefly, approximately 10–15 g of sieved (< 2 mm), air-dried soil was placed in weighed Büchner funnel with pre-wetted Whatman No. 2 filter paper, and then, the soil was wetted repeatedly with deionized water and allowed to drain by gravity for 48 h at 100% humidity in a closed environment. Vacuum suction (− 33 kPa) was then applied to drain the soil to approximate field capacity prior to reweighing. We calculated soil WHC as the water mass retained divided by oven-dry soil mass equivalent.

Particle size distribution (PSD; % sand, silt, and clay of the < 2 mm soil fraction) was determined using the hydrometer method (Bouyoucos 1962; Gee and Bauder 1986). For PSD analysis, we repeatedly added 30% hydrogen peroxide as a pre-treatment for all soil samples to remove organic matter particles, as several meadow soils had relatively high organic matter content. Soil organic matter (SOM) content was determined by loss on ignition (modified from Combs and Nathan 1998) from sieved (< 2 mm) soil subsamples placed in a muffle furnace (360 °C for 2 h).

Laboratory analyses were conducted in two groups (group A and B) to reduce laboratory time and cost. Some analyses (BD, GWC, RMAD, percent roots, and coarse fragments) were done for all soil cores (group A, n = 252), while other analyses (SOM content, PSD, WHC) were done on only a subset of soil cores (group B, n = 98).

Statistical analysis

After laboratory analyses were conducted on the two soil core groups (group A and B), the groups were stratified by either plant community type (PCT) or meadow-scale gradient class (MGC) in order to separate the samples by spatial scale. Equal sample sizes were ensured among strata for both groups, as shown in the stratified random study design diagram schematic (Fig. 2). RStudio software (R Core Team 2017) was used to conduct statistical analyses. For all statistical tests conducted, statistical significance was set a priori at α = 0.05.

To address our research questions, we examined the covariance of potential explanatory factors with SR at the local scale (among PCTs) and at the meadow scale (between MGCs). We conducted this assessment to test whether vulnerability to trampling by pack stock animals in localized plant communities within meadows is contingent on the meadow scale (among meadows of differing topographic gradients).

To address our first research question, we used Pearson’s product moment correlation analysis (“corrplot” package in R) to identify potential explanatory factors that best correlated with SR and to identify and avoid combining collinear explanatory variables (r ≥ 0.75) in our linear mixed-effects regression (LMER) models (Graham 2003). Multiple LMER models (“lme4” package in R) were used to determine which combination of potential explanatory factors best explained variation in SR. The LMER mixed modeling approach accounts for a stratified study design (Zuur et al. 2007) and takes into account fixed variability among PCTs and meadow gradient classes (MGCs) and random variability among the study meadows (Bates et al. 2012; Roche et al. 2014). Our use of mixed effects modeling assumes that our study meadows are a random subset of a large and variable population (Gruebber et al. 2011; Holmquist et al. 2014; Roche et al. 2014), which allows greater confidence when inferring our results to other meadows not sampled in our study. Our study included a subset of five meadows within a total population of over 3000 meadows in Yosemite, which have variable plant and soil characteristics. We used Akaike and Bayesian information criterions (AIC and BIC) to identify the best explanatory model, balancing goodness-of-fit with parsimony, from a range of potential models (Crawley 2007; Bates et al. 2012). Our null models were comprised of only PCT or MGC grouping factors or both, where “meadow” was identified as a random factor (for each of the five meadows). We identified best models as those with AIC values lower than 10 points or more than the null model and with the lowest corresponding BIC values (Dziak et al. 2012; Burnham and Anderson 2002). We also used diagnostic tests to check model assumptions (e.g., normality, linearity, independence, and non-constant variance); no additional transformations were needed to meet these model assumptions.

To address our second research question, we used a three-factor analysis of variance (ANOVA) model to represent the two spatial scales of our study: local-scale plant community types (PCTs, n = 4) and meadow-scale gradient classes (MGCs, n = 2), and temporal variation among study weeks (n = 3). Recall that soils were sampled in plots over 3 weeks due to the remoteness of study locations, where only a small number of samples could be carried out of the field weekly. Hence, we chose to also evaluate temporal variation among study weeks to test whether timing of sampling may have affected our results, although the main focus of our study is on spatial variation at two scales. We used model interaction terms to evaluate whether the patterns of SR and explanatory variables are contingent on the spatial scale (interaction between PCT and MGC) or timing of measurement (interaction with week). We assessed normality using histograms (“car” package in R) and tested homogeneity of variance (HOV) using HOV plots (“HH” package in R). To achieve normality and homogeneity of variance, we natural log-transformed the SR response variable and square-root transformed the GWC, RMAD, root content, and SOM content potential explanatory variables; no other explanatory variables required transformation.

On a local scale, we also tested whether PCTs differed in the ability to support pack stock use without soil deformation by comparing mean SR values by to a threshold of 500 kPa using one-way ANOVA (n = 252). This threshold value is a conservative estimate that represents the SR needed to support a horse with rider or mule with load without incurring soil compaction, based on modifications to (Scholefield and Hall 1985) and Kai et al. (2000). We modified these approaches, by tailoring our kPa threshold to that of pack stock, rather than livestock animals.