Multi-scale drivers of soil resistance predict vulnerability of seasonally wet meadows to trampling by pack stock animals in the Sierra Nevada, USA

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. 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). 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. 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.


Introduction
Meadows in the western United States of America provide ecosystem services and socioeconomic benefits, such as high plant diversity, critical habitat for wildlife, sediment and water storage and filtration, nutrient cycling, and flood attenuation, and are also often popular recreational destinations Russo et al. 2012;Norton et al. 2014). Meadows contribute these services in far greater proportion than their very limited areal extent of less than 1% of the Sierra Nevada and Cascade Range (Viers et al. 2013) and only 3% of Yosemite National Park (Keeler-Wolf et al. 2012). On public lands, high elevation meadows are popular destinations for day hikers, backpackers, and recreational or administrative users with pack stock animals (horses and mules) due to their scenic beauty, iconic mountain vistas, close proximity to water, and availability of high-quality summer forage. Yet, these important ecosystems are vulnerable to degradation, especially when soils are wet immediately after snowmelt and vegetation is first developing in early summer (Cole et al. 1987). In the Sierra Nevada, the potential for seasonally wet meadow soils to incur damage from these human uses is a main concern for land managers (DeBenedetti and Parsons 1983;Cole et al. 1987;Shryock 2010;Lee et al. 2017;Ostoja et al. 2014).
Soil resistance (SR), as a measure of soil strength, has been used to quantify the susceptibility of soil to physical disturbance by humans and other animals (Lull 1959;Herrick and Jones 2002;Herbin et al. 2011). Soil resistance is defined as the ability of soil to resist deformation under applied pressure; this soil property varies inversely with soil water content because disaggregation of soil particles increases with increasing soil water content (Scholefield and Hall 1985;Perumpral 1987;Herrick and Jones 2002;Piwowarczk and Holden 2011). Conversely, SR varies directly with soil bulk density because soil strength increases as inter-particle forces and friction increases (Perumpral 1987;Herrick and Jones 2002;Herbin et al. 2011). Increases in SR occur with increases in bulk density (i.e., soil compaction) and associated reductions in soil water holding capacity (Herrick and Jones 2002;Piwowarczk and Holden 2011).
Many ecological factors influence the soil water content in seasonally wet mountain meadows, where snow water equivalent and spring snowmelt timing influence duration of high-water tables, length of soil saturation, timing of plant phenology, and plant productivity (Loheide et al. 2008;Loheide and Lundquist 2009;Moore et al. 2013). Plant community composition is strongly correlated to soil water content at both the local scale (within a meadow) and meadow scale (the overall meadow; Ratliff 1985). At the local scale, soil water content (and also soil redox potential) is influenced by soil structure, water table depth and dynamics, and soil organic matter content (Chambers et al. 1999;Castelli et al. 2000;Rodríguez-Iturbe and Porporato 2004;Dwire et al. 2006;Loheide and Gorelick 2007;Loheide et al. 2008). At the meadow scale, soil water content is influenced by hydrogeomorphic type and topographic gradient (Brady and Weil 2008;Huber et al. 1989;Weixelman et al. 2011;Roche et al. 2014). The hydrogeomorphic approach for meadow-scale classification accounts for functional differences among diverse meadow and wetland types; meadows are classified by several characteristics, such as the presence of a stream channel (riparian) and steepness of the topographic gradient (high, middle, or low), which is determined by the overall meadow slope (Brinson 1993;Weixelman et al. 2011).
Our study aimed to determine whether multi-scale drivers of soil resistance (SR) can be used to predict wet meadow vulnerability to trampling by pack stock animals in the Sierra Nevada. Specifically, we investigated which factors best explained local-scale variation in SR and whether local-scale variation in drivers of SR was contingent on meadow-scale gradient class. We used SR as an index of the vulnerability of undisturbed soils to compressional deformation (trampling) by pack stock animals. Our study addressed the following two research questions: (1) What factors best explain localscale variation in SR (among plant community types in meadows)? and (2) Are local-scale (plant community-level) differences in SR and potential explanatory factors of SR contingent on the meadow-scale gradient classes? To address our study questions, we measured SR and several soil and vegetation properties at two spatial scales using a stratifiedrandom sampling design (in dominant plant community types within seasonally wet subalpine meadows representing low and middle topographic gradient classes). Plant community classification was used to recognize within-meadow (local scale) variation (Ratliff 1985;Allen-Diaz 1991;Loheide and Gorelick 2007;Loheide et al. 2008), and meadow hydrogeomorphic and gradient classification was used to recognize among-meadow (meadow scale) variation (Huber et al. 1989;Brady and Weil 2008;Weixelman et al. 2011;Roche et al. 2014). At the local scale, we expected that SR would be related primarily to soil texture and resulting soil water content. At the meadow scale, we expected that variability in SR would be most strongly influenced by gradient class and resulting soil texture and water content. Hence, our overall hypothesis was that local-scale variation in drivers of SR would be contingent on meadow-scale gradient class, indicating that multi-scale drivers of SR can be used to predict wet mountain meadow vulnerability to trampling by pack stock animals.

Study area
Our study was conducted in Yosemite National Park (Yosemite) in the central Sierra Nevada of California, within subalpine meadows in the Tuolumne watershed. This region experiences a Mediterranean-type climate, with warm, dry summers and cool, moist winters ( Fig. 1). Tuolumne Meadows, one of the study sites representative of all the meadows studied, has a mean annual precipitation of 755 mm, of which 80-90% falls as snow generally between October and April. Mean daily air temperatures vary from − 12.7°C in January to 21.3°C in July, and snowmelt occurs in a large pulse typically in May and June (Clow et al. 2010;Moore et al. 2013). Our study was conducted during peak growing season (July) in one of the driest years on record in the Sierra (2013). The 2013 average April 1 snow water equivalent (SWE) was 50% of average, which was the lowest SWE during a prolonged drought period (DWR, CDEC 2019).
This study was part of an interdisciplinary meadow investigation conducted in Yosemite on seasonally wet meadow vulnerability to stock use impacts, which investigated a total of 65 meadows subject to pack stock use in the Tuolumne watershed (Kuhn et al. in review). In this study, we did not evaluate the effects of grazing on meadow soils, nor did we evaluate grazing intensity or duration. Please see (Jones et al. 2018) for information on meadow grazing capacities of Yosemite National Park. Of the~3000 meadows within Yosemite Wilderness (Viers et al. 2013), only 200 are within the 0.4 km (1/4 mile) limit prescribed by park regulations and policies for pack stock use, where only 65 meadows had recorded grazing since 2004 (Kuhn et al. in review). Within the study area, we selected five meadows within the subalpine zone (2450 to 3250 m) that are located at similar elevations (2667 to 2849 m) and with similar soils (i.e., Xeric and Oxyaquic Dystrocryepts; NRCS 2006). These meadows represented two different meadow gradient classes according to the hydrogeomorphic approach for meadow-scale classification (Weixelman et al. Points on map depict study meadows used to evaluate explanatory factors of soil resistance by hydrogeomorphic type and appearance (green circles: low gradient meadows, < 2% slope, yellow squares: middle gradient meadows, 2-4% slope) 2011): low gradient, slopes < 2%; and middle gradient, slopes 2-4%. We calculated mean percent meadow slope (meadow scale) using a 10-m resolution (horizontal precision) digital elevation topographic model. Low gradient meadows included Snow Flat (1.5%, 4 ha), Tuolumne Meadows-West (1.8%, 38 ha), and Upper Lyell Canyon (1.9%, 14 ha), and middle gradient meadows included Middle Lyell Canyon (3.7%, 6 ha) and Emeric Lake (3.0%, 10 ha). All meadows are riparian except Emeric Lake, which is partially lacustrine fringe (i.e., near a lake). All meadows are also subject to pack stock use, except Tuolumne Meadows-West, which represents the largest meadow in the park and is a popular day-use destination.
We selected four plant community types (PCT), representative of Sierran meadows, based on botanical study data for monitoring and management of pack stock use in Yosemite high elevation wilderness meadows (Kuhn et al. 2015). These PCTs are indicative of different hydrological meadow types that represent a soil moisture gradient from wet to dry within meadows (Baldwin et al. 2012;Keeler-Wolf et al. 2012). Specifically, we assigned soil moisture types, or hydrologic regimes (xeric, mesic, hydric), to the plant community types (PCTs) we studied based on wetland indicator status, as defined by USACE (2013): hydric PCT-inflated-beaked sedge (Carex vesicaria; CV), based on an obligate (OBL) wetland status; Mesic PCT 1-tufted hairgrass (Deschampsia cespitosa; DC), based on a facultative-wetland (FACW) status; Mesic PCT 2-Muir's reedgrass (Calamagrostis muiriana; CM), based on a facultative-wetland (FACW) status; and Xeric PCT-King's ricegrass (Stipa kingii; SK), based on a facultative upland (FACU) status.

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 2 plots placed in four . Soil and vegetative cover sampling occurred over three study weeks in July of 2013: July 15th (sample 1), July 22nd (sample 2), and July 29th (sample 3). Weekly sampling occurred within randomized plots within PCTs. *Plot sample analysis was done for the soil resistance response variable and 14 explanatory variables: group A (n = 252, derived from sampling six plots with two combined replicates per plot, per PCT)-soil resistance, soil bulk density, soil gravimetric water content, soil volumetric water content, soil coarse fragments 2-4 mm in diameter, soil coarse fragments > 4 mm in diameter, total vegetation cover (at 1 m 2 and 700 cm 2 scales), root mass areal density, and root content; group B (n = 98)-soil water holding capacity, soil organic matter content, sand content, silt content, and clay content meadow PCTs (i.e., strata; Fig. 2). 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 (Fig. 2).
We used meadow gradient class to examine meadowscale 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 2 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 2 sample plots, we used a nested sample-plot design of 15-cm and 5-cm radius circular plots (~700 cm 2 and~80 cm 2 , respectively) to measure SR, soil volumetric water content (VWC), and total vegetation cover (TVC). In the 1-m 2 and~700-cm 2 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~80cm 2 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 crosssectional area of the core (19.6 cm 2 ).
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.

Results
Factors that best explain local-scale variation in soil resistance among meadow plant community types Eight of 14 explanatory variables were significantly correlated with SR (r > ± 0.16, p < 0.05 for n = 252 and r > ± 0.20, p <0.05 for n = 98; Table 1). The strongest significant relationship was a negative correlation with GWC (r = − 0.33), closely followed by a positive correlation with BD (r = 0.31). Three sets of explanatory variables were determined to be collinear (r > 0.75, Table 1): GWC was collinear with WHC (r = 0.88) and BD (r = − 0.86), and RMAD were collinear with root content (r = 0.87). To avoid inflation of model explanatory power (Graham 2003), we did not combine collinear variables in any one model when determining best explanatory models of SR.
Using LMER analysis, we identified five best explanatory mixed models based on non-collinear, parsimonious combinations of explanatory variables with the lowest AIC and BIC model criterion (Table 2). We compared best mixed models with associated combinations of null and single-variable models. Null models included one or more grouping factors (PCT, MGG) and a random factor only (Mdw). Single-variable models included the random factor (Mdw) and one main explanatory factor only (i.e., BD, GWC, or RMAD). Best mixed models included both grouping factors (PCT and MGC), the random factor (Mdw), and a combination of up to two of the following three explanatory variables: BD, GWC, and RMAD. The best explanatory models were mixed model 1 (BD and RMAD, main factors, AIC = 209, BIC = 232), mixed model 2 (BD main factor, AIC = 211, BIC = 232), mixed model 3 (GWC and RMAD main factors, AIC = 210, BIC = 233), and mixed model 4 (GWC main factor, AIC 210, BIC = 230). The best mixed model (mixed model 1) had an AIC value that was at least 16 points lower than null model AIC values (AIC range = 225-244) and at least 10 points lower than null model BIC values (BIC range = 243-254). When evaluating both AIC and BIC together, the most parsimonious model was a single-variable model (GWC main factor), which had the second lowest AIC score (AIC = 210) and the lowest BIC score (BIC = 220). This single-variable model had an AIC value that was at least 15 points lower than the null model AIC values (AIC range = 225-244) and at least 23 points lower than null model BIC values (BIC range = 243-254).

Local-scale variation in drivers of soil resistance is contingent on meadow-scale gradient class
We sampled up to three of the four PCTs selected for study in each of our five study meadows because all four plant community types do not typically occur together in any one meadow. We found that local-scale Table 1 Pearson's product-moment correlation (r) matrix results for the soil resistance (SR) response variable and 14 explanatory variables. Explanatory variables derived from soils sampled over a 3-week period in subalpine meadows (n = 5) at Yosemite National Park, USA: group A (n = 252) bulk density (BD), gravimetric water content (GWC), volumetric water content (VWC), total vegetation cover (TVC 1 m 2 ), TVC 700 cm 2 , coarse fragment content (CF) 2-4 mm, CF > 4 mm, root mass areal density (RMAD), and root content (roots %), group B (n = 98): water holding capacity (WHC), soil organic matter content (SOM), and sand, silt, and clay. Bolded values represent variables significantly correlated with SR (r > ± 0.16, p < 0.05 for n = 252 and r > ± 0.20, p < 0.05 for n = 98), or significantly correlated with each other, indicating collinearity (r ≥ ± 0.  (Tables 3  and 4; Fig. 3). During the study period, neither of the wet plant community types (DC nor CV) had SR values that exceeded the threshold needed to support a horse or mule (i.e., 500 kPa, Fig. 3). The mean values for the response variable and each explanatory factor are provided in the supplementary materials (Tables S1-S9). Several potential explanatory variables also showed significant interactions between PCT and MGC. For instance, BD was generally lower in the wetter PCTs as follows: CV<SK<CM and CV<DC<SK in low gradient meadows, but CV<DC<CM<SK in middle gradient meadows. Soil VWC was generally higher in the wetter PCTs, contingent on MGC as follows: SK<CM<DC and SK<CM<CV in low gradient meadows, but SK<CM< DC<CV in middle gradient meadows. Small coarse fragments (CF 2-4 mm) were generally greater in the dry PCTs and middle gradient meadows. Total vegetative cover in the 1-m 2 plots ranked CV<SK<CM and CV< SK<DC in the low gradient meadows, but ranked CV< CM<DC and CM<SK in the middle gradient meadows. Roots (%) also exhibited a significant PCT by MGC interaction; roots (%) ranked CM<SK<DC and DC<SK< CV in low gradient meadows, while roots (%) ranked CM<DC<CV and CV<SK in middle gradient meadows.
Large coarse fragments (CF > 4 mm) also exhibited a significant PCT by MGC interaction, but this interaction depended on the week (which covaried with the plot) sampled.
Only four variables exhibited no statistically significant two-or three-factor interactions in the three-factor ANOVA models (Table 3). Gravimetric water content (GWC) was significantly affected by PCT (p = 0.00) and MGC (p = 0.00) factors. Sand, silt, and clay content (all soil texture related) were significantly affected by one or both grouping factors. Sand and clay were affected by both PCT (p = 0.03 and 0.02, respectively) and MGC (p = 0.00 for both). Clay content was significantly affected by MGC (p = 0.00) but not PCT (p = 0.44). In low gradient meadows, patterns in GWC showed higher soil water content in wetter PCTs: CV>DC>CM>SK. In middle gradient meadows, a slightly different pattern in GWC was detected among PCTs: CV>DC>SK>CM. Patterns in % sand and % silt were inverse of each other, where low gradient meadows had more silt content and less sand content and middle gradient meadows had more sand content and less silt content. Sand and silt contents among PCTs in low gradient meadows were as follows: PK>DC>CM>CV and CV>DC>CM>PK, respectively. In middle gradient meadows, sand and silt content did not show the same patterns: PK>DC>CV>CM and CM>CV>DC>PK, respectively. Overall patterns in GWC, % sand, % silt, and % clay among PCTs indicate that these variables are not Table 2 Linear mixed-effects regression (LMER) results showing best explanatory models of soil resistance (SR) in comparison with null models. Best models are shown in bold font, based on lowest Akaike (AIC) and/or Bayesian (BIC) Information Criterions. Null models include one or more grouping factors and a random factor only, single models include one main factor (explanatory variable) and random factor only, and best models include one or more grouping factor, the random factor, and one or more main factors for dataset group A (n = 252) derived from soils sampled over 3-week period in in subalpine meadows (n = 5) at Yosemite National Park, USA. An "X" denotes model factors included within a given model. Grouping factors included plant community type (PCT) and meadow gradient class (MGC), and the random factor is each of the five subalpine meadows (Mdw). Main factors included the following explanatory variables: bulk density (BD), gravimetric water content (GWC), and root mass areal density (RMAD)

Models
Grouping contingent on meadow gradient (i.e., no interaction occurred between spatial scales; Tables 3 and 4). While GWC, % sand, and % silt significantly differed at both spatial scales, % clay significantly differed only at the meadow scale (Tables 3 and 4; Tables S1-S9).

Discussion
At the local scale within meadows, we found that SR was related primarily to soil texture and resulting soil water content. Three explanatory variables, representing soil disaggregation, friction, and cohesion, best explained local-scale variation in soil resistance (SR) among meadow plant community types. Based on our explanatory and null models using linear mixed effects regression (LMER), these variables were gravimetric water content (GWC), bulk density (BD), and root mass areal density (RMAD). Specifically, GWC covaries with soil disaggregation, and BD and RMAD covary with soil inter-particle friction and cohesion (Brady and Weil 2008). Previous livestock-use studies in grassland ecosystems have shown that soil water content varies inversely with SR because disaggregation of soil particles increases with increasing soil water content (Sholefield and Hall 1985;Perumpral 1987;Herrick and Jones 2002;Piwowarczk and Holden 2011). These earlier studies also found that SR varies directly with BD because soil strength increases as inter-particle forces and friction increases. They concluded that both GWC and BD best explain variation in SR (Perumpral 1987;Herrick and Jones 2002;Herbin et al. 2011). Our results support and also build upon these previous investigations, showing that RMAD is positively correlated to local-scale SR among meadow PCTs, likely due to roots' positive impact on soil friction and cohesion (Brady and Weil 2008). We found that the most parsimonious explanatory model of SR (based on lowest AIC and BIC criteria) included GWC alone, closely followed by the model including BD alone. The mixed model that included BD and RMAD as main factors was the best model based on AIC score alone, but only by a factor of one point. The improvement for using RMAD as a second variable Table 3 Results of three-factor ANOVAs for soil resistance (SR) and various explanatory variables. Two grouping factors represent different spatial scales: plant community type (local scale) and meadow gradient class (meadow scale). The other grouping factor (week) denotes the sampling week that the variable was measured over a 3-week study period (temporal effect). Note: different plots were sampled each week across all five meadows. Bolded F and p values (F) indicate variables with significant differences (p < 0.05) for that factor or interaction. Acronyms and abbreviations: BD bulk density, GWC gravimetric water content, VWC volumetric water content, WHC water holding capacity, SOM soil organic matter content, Sand % sand content, Silt % silt content, Clay % clay content, CF coarse fragment content in two size classes (2-4 mm and > 4 mm), TVC total vegetation cover in two plot sizes (1 m 2 and 700 cm 2 ), RMAD root mass areal density, Roots % root content differed depending on the initial variable used; it improved the AIC score for the BD-only model, but did not improve the AIC score of GWC-only model. We hypothesize that BD and RMAD have complimentary influence on aggregation and soil structure, and thus soil strength. This complimentary influence may be the result of combining variables with different spatial scales of influence: BD differed among local-scale PCTs, but SR and RMAD among PCTs were contingent on the meadow scale (MGC). Our results suggest that SR is determined primarily by soil texture and resulting water content at the local scale within meadows and that meadow-scale variability in SR is most strongly influenced by topographic gradient and resulting soil texture and water content. Low gradient meadows surveyed had finer substrates, greater GWC, wetter PCTs, and lower SR, resulting in greater soil disaggregation and higher vulnerability to soil disturbance. Conversely, middle gradient meadows surveyed had coarser substrates, lower GWC, drier PCTs, and higher SR, resulting in greater soil friction and cohesion and lower vulnerability to disturbance. Hence, our data suggest that Sierran meadows with middle gradients (2-4% slope) have lower vulnerability to trampling damage by pack stock use than lower gradient meadows (< 2% slope). Similarly, we speculate that high gradient meadows (> 4% slope), not assessed in this study, would also have low vulnerability to damage, based on our findings for middle gradient meadows.
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 (Holmquist et al. 2014). Our study found that spatial patterns of SR are dependent on the scale of observation, where local-scale differences among PCTs were different between the two meadow-scale gradient classes (MGC) studied. The positive correlation between SR and bulk density (BD) indicates that local-scale differences in soil strength (among PCTs) were also contingent on the meadow scale (between MGCs), in that both of these variables had a significant interaction between PCT and MGC grouping variables (Table 3). At the meadow scale, we found that variability in SR was most strongly influenced by topographic gradient and resulting soil texture and water content. Earlier livestock studies in grasslands also reported consistent local-scale spatial patterns between SR, BD, and water content (Perumpral 1987;Herrick and Jones 2002;Herbin et al. 2011), which suggests similar controls on SR across herbaceous ecosystems. In addition to soil BD and water content, we found significant relationships between SR and several other properties in meadow ecosystems, including vegetation cover, root mass, and coarse fragments, as well as a significant relationship with MGC. For example, we found that large coarse fragment content corresponded with greater SR and lower VWC. This is consistent with the findings of Childs and Flint (1990) who found that increases in rock fragment content in soil can cause dramatic changes in ecosystem processes and properties, including decreases in plant productivity and decreases in soil silt and clay content, water content, and nutrient content. These overall findings provide a better understanding of how meadow vulnerability to pack stock trampling varies on the local and meadow scales.
We found that neither of the wet plant community types (DC nor CV) had SR values that exceeded the threshold (500 kPa) needed to support a horse or mule during the study period. This was the case regardless of the MGC over the entire study period. Exposure to pack stock use in wet PCTs may seem unlikely in wet areas; however, a recent Yosemite stock behavior study showed that horses and mules graze on Carex vesicaria (the wettest PCT in our study), which makes this plant community type vulnerable to trampling when soil moisture conditions are high (Walden-Schreiner et al. 2017). In contrast, we found that local-scale plant community vulnerability to soil physical damage was modulated by Table 4 Spatial interactions between plant community type (PCT, n = 4) and meadow gradient class (MGC, n = 2) within subalpine meadows sampled (n = 5) that were significant based on ANOVA (n = 252 for all variables). Rows display relative ranking of values for that variable for each PCT that were contingent on MGC. Plant abbreviations and wetland status: CV Carex vesiciaria (obligate), DC Deschampsia cespitosa (facultativewetland), CM Calamagrostis muiriana (facultative), SK Stipa kingii (facultative-upland). Acronyms and abbreviations: SR soil resistance, BD bulk density, VWC volumetric water content, CF coarse fragment content in two size classes (2-4 mm and > 4 mm), TVC total vegetation cover determined in 1 m 2 plots (1 m 2 ), RMAD root mass areal density, Roots % root content meadow-scale gradient class in drier PCTs. For instance, in the xeric PCT (SK), middle gradient meadows had sufficient SR to support pack stock without soil deformation, but this was not the case for this xeric PCT in low gradient meadows (Fig. 3). Due to the lower surface water kinetic energy across the landscape (i.e., sheet flow hydrology), low gradient meadows have finer soil textures, higher soil organic matter contents, and higher water holding capacities (Brady and Weil 2008;Weixelman et al. 2011;Viers et al. 2013;Wang et al. 2013;Roche et al. 2014), resulting in higher vulnerability to trampling.
Our study sampled only five subalpine meadows out of approximately 18,780 meadows in the Sierra Nevada (UCD SNMDC 2019). Additionally, due to logistical constraints related to study site remoteness, we sampled meadows over a 3-week period in order to incorporate enough replicate plots to provide statistical inference. Climate change may also be an important factor to consider when evaluating the patterns we observed. For instance, our study was conducted during peak growing season (July) in one of the driest years on record in the Sierra (2013). The most likely impact of drought on our study is that the time needed to achieve sufficient soil strength to resist soil deformation from pack stock impacts was much shorter than it would be in a year with average or above-average precipitation. Importantly, we found that even during drought, SR in wet plant communities was insufficient to sustain pack stock use without incurring soil deformation. Hence, we speculate that wetter years would result in even greater delays in achieving sufficient soil strength to minimize pack stock trampling impacts in these meadows. In wetter years, low gradient meadows with wet plant communities may never achieve sufficient soil strength needed to resist soil deformation from pack stock use.

Conclusions
We confirmed our overall hypothesis that local-scale variation in drivers of soil resistance (SR) is contingent on meadow-scale gradient class, indicating that multiscale drivers of SR can be used to predict wet meadow vulnerability to trampling by pack stock animals. Our study provides new information on the multi-scale drivers of SR within high-elevation meadows and should be useful for developing a predictive tool for assessing Fig. 3 Mean soil resistance (SR) by plant community type (PCT) among meadows sampled over a 3-week study period in July, 2013 (n = 252). Dashed line indicates SR minimum value (500 kPa) needed to support a horse with rider or mule with load without incurring physical soil disturbance. Error bars represent plus one standard error (SE) across plots within a PCT and meadow. X-axis displays meadows of two topographic gradient classes. Middle gradient meadows (left): Emeric Lake (EL) and Middle Lyell (ML). Low gradient meadows (right): Snow Flat (SF) Tuolumne Meadows (TM), and Upper Lyell (UL). Missing bars indicate a PCT was not present in that meadow seasonally wet meadow vulnerability to trampling from pack stock use in the Sierra Nevada and elsewhere where meadows are dominated by similar herbaceous plants with seasonal dry-down patterns. We build upon earlier studies highlighting the importance of soil water content and bulk density values in determining soil resistance to deformation by pack stock use. In addition to these controls, we also found that root mass areal density, apparently through its influence on soil friction and cohesion, also influences SR. Our study advances the development of a meadow vulnerability predictive tool that can help inform meadow management of pack stock in seasonally wet, high-elevation meadows in the Sierra Nevada. While small in extent, wet meadows represent some of the last remaining intact wetlands in the state, and they store the source of a large amount of the state's water supply (snowmelt), making them critically important to the socioeconomic services they support (Viers et al. 2013). Often located on public lands, high-elevation meadows are increasingly at risk to recreational use, making it imperative that these ecosystems be managed wisely with effective ecosystem-based approaches (Ostoja et al. 2014). Sierra Nevada meadows are highly degraded throughout the range (Viers et al. 2013), and our study should be useful to resource managers for making scientifically sound management decisions on meadow use by pack stock animals.