Elevation gradient drives distribution of soil carbon in a semiarid grassland of British Columbia

Abstract A sequence of Brown, Dark Brown and Black Chernozems spanning a 600 m elevation gradient in a semiarid bunchgrass ecosystem (Lac du Bois Grassland) near Kamloops, British Columbia was first described in 1961. More soil organic carbon (SOC) at higher elevations along the sequence was attributed to increasing effective precipitation with increasing elevation. Since the 1961 study, plant community composition has shifted toward the desired climax community due to improved livestock management instituted in the 1970s; however, changes in soil carbon stocks remain unknown. The objective of this study was to quantify SOC and soil inorganic carbon (SIC) stocks using the same site selection criteria as used in 1961. SOC stocks (kg m−2 ± SD; 0–60 cm) were similar for Brown (5.73 ± 1.7) and Dark Brown Chernozems (5.87 ± 0.76) but increased sharply (10.11 ± 2.5) for the higher elevation Black Chernozems. SIC increased with depth in all three soil zones, representing 33%–50% of total C from the 30–60 cm soil depth. To evaluate changes in SOC (0–20 cm) from the 1961 measurements, three different approaches for calculating SOC stocks were used based on the inclusion or exclusion of coarse fragments. Results varied across the three soil zones from no change to a 20% increase in the Brown, an increase of 7% to a reduction of 26% in the Dark Brown, and a decrease of 12% to 35% in the Black soil zone. Information about soil coarse fragments and the distribution of SOC and SIC stocks within the soil profile is crucial for accurate comparisons across studies or resampling events.


Introduction
Grassland ecosystems cover more than 40% of the global terrestrial area (Hewins et al. 2018;Bardgett et al. 2021) and offer numerous ecosystem services such as wildlife habitat, livestock forage, and soil organic carbon (SOC) storage, which is estimated to represent up to 30% of the global terrestrial C pool (Abdalla et al. 2018;Hewins et al. 2018;Bai and Cotrufo 2022).Grassland soils are a critical C sink in the global C cycle due to the ongoing addition and stabilization of organic matter that typically outweighs the slow decomposition due to limited soil water (Pennock et al. 2011).However, historical and present-day grazing management have degraded up to 30% of grasslands worldwide, representing a loss of SOC between 1.8 and 6.3 Gt C (Dlamini et al. 2016).This emphasizes the need to implement grazing strategies that increase plant cover and encourage plant species that maximize SOC storage.The mechanisms that control SOC sequestration in grassland soils differ by climate (Conant et al. 2017;Sanderson et al. 2020), soil type (Dlamini et al. 2016;Reinhart et al. 2021), grazing intensity (Han et al. 2008), and photosynthetic pathway of the dominant grassland species (McSherry and Ritchie 2013;Abdalla et al. 2018;Liang et al. 2021).A better understanding of soil C stocks at a regional level will greatly improve our ability to better predict C storage capability of grasslands on a global scale (Hewins et al. 2018).
Organic matter in grassland soils is primarily mineralassociated organic matter (MAOM; <50 μm) with long turnover times (Cotrufo et al. 2019;Lavallee et al. 2020).Leaf litter inputs contribute mostly to particulate organic matter (POM) with shorter turnover times (Steffens et al. 2011;Li et al. 2022), whereas root inputs contribute more to MAOM formation and longer C sequestration (Bai and Cotrufo 2022).Therefore, grazing regimes that support plant species that allocate more C belowground will favour greater SOC storage.However, grazing-induced changes on SOC may be obscured by the heterogeneous nature of organic matter compounds, varying degrees of plant and root litter recalcitrance, and chemical and physical protection of those compounds with soil mineral particles (Schmidt et al. 2011).In a study from Inner Mongolia, MAOM fractions did not increase after 25 years of grazing cessation in semiarid grassland because soil particles <50 μm were thought to be saturated and the increased aboveground biomass returned to the soil was primarily contributing to POM fractions, rather than MAOM (Steffens et al. 2011).As a result, POM is considered to be a more sensitive indicator of grazing-related change in SOC (Dong et al. 2021).Grassland soils are also known to have high concentrations of carbonates (e.g., CaCO 3 ) that contribute to a substantial soil inorganic carbon (SIC) pool (Heckman and Rasmussen 2018).Because total C measurements will include both the organic and inorganic C pools, the amount and location of pedogenic carbonates (i.e., CaCO 3 ) within the soil profile should be explicitly measured so that SOC can be accurately accounted.Furthermore, soil texture has been a strong predictor of both SOC and SIC, due to the influence that parent material has on the distribution of carbonates (i.e., SIC) and the ability of fine soil particles (<50 μm) to absorb more soil C (Deane McKenna et al. 2022).Therefore, the potential of a particular soil to sequester additional soil C becomes a function of soil texture (Georgiou et al. 2022), the formation and persistence of POM (Six and Paustian 2014), and the degree of MAOM saturation (Stewart et al. 2009;Georgiou et al. 2022).
The majority of Canadian grasslands are found in the Prairie provinces, with only a small portion in British Columbia (BC) (Pennock et al. 2011).Grassland ecosystems and associated Chernozems in BC are located mainly in the Central and Southern Interior, the Peace River region, Cariboo-Chilcotin, and Rocky Mountain Trench bordered by the Rocky Mountains to the east and Columbia Mountains to the west (Wikeem and Wikeem 2004).Most of BC's grasslands are associated with an undulating topography with significant elevation changes over short distances, in turn leading to a diversity of microclimates, plant communities, and soil types.One of the most studied grasslands of BC is the Lac du Bois Grassland, located just north of Kamloops (Supplemental S1).
Soil and vegetation of the Lac du Bois Grassland were initially described by Spilsbury and Tisdale (1944) through a series of soil-climatic zones over a 10 km distance that change along an elevation gradient (400-1000 m above sea level (asl)).The corresponding increase in SOC across the sequence of Brown, Dark Brown, and Black Chernozems on south-facing mid-slope positions was described in 1961 by van Ryswyk et al. (1966) as a function of increasing precipitation and decreasing temperature with elevation gain.Like most grasslands in BC, the Lac du Bois Grassland has been severely degraded by heavy grazing during late 19th to mid-20th centuries.With the introduction of improved grazing strategies throughout south-central BC since the mid-1970s, that either reduce the grazing duration or alternate the season of grazing, there has been an increase in perennial bunchgrass cover and minimal differences in species richness and diversity between grazed and ungrazed sites (Bradfield et al. 2021).Specifically at Lac du Bois, livestock grazing is managed using a seasonal deferred grazing regime (Briske et al. 2008).Beginning in April, livestock are introduced to alternating pastures in every year and grazed until June.No grazing occurs during summer (June-September), but grazing returns to the grassland in October and continues until snowfall.Although improved grassland condition with better management has been well known, comprehensive assessments of the SOC stocks are still lacking.The objective of this study was to quantify SOC and SIC stocks (kg m −2 ) along the elevation gradient using the same sample site se-lection criteria used by van Ryswyk et al. (1966) since the introduction of an improved grazing strategy at the Lac du Bois Grassland.

Study site
The Lac du Bois Grassland, located just north of the city of Kamloops, BC (latitude 50.757, longitude −120.418), is classified as a dry-steppe ecosystem (Carlyle et al. 2011) within the interior dry belt of BC and the Montane Cordillera Ecozone.Daubenmire (1970) referred to this area as the most northern reaches of the Palouse Prairie with close associations with the Pacific Northwest Bunchgrass ecosystem that extends through portions of Oregon, Washington, Idaho, and Montana (Endress et al. 2012).Three plant communities were described at Lac du Bois in terms of their elevation sequence by Spilsbury and Tisdale (1944): (i) Agropyron-Artemisia or Lower Grassland Zone, (ii) Agropyron-Poa or Middle Grassland Zone, and (iii) Agropyron-Festuca or Upper Grassland Zone.The interaction of temperature and precipitation effectiveness best described the relatively sharp boundaries between the grassland zones (van Ryswyk et al. 1966).Using a similar distinction between zones, the Lower, Middle, and Upper grasslands can also be classified as the Bunchgrass Very Dry Hot Thompson variant (BGxh 2 ), Very Dry Warm Nicola variant (BGxw 1 ), and as a grassland phase within the Interior Douglas-fir Very Dry Hot Thompson variant (IDFxh 2 ), respectively (Ryan et al. 2021).Dominant plant species in the Lower Grassland Zone are bluebunch wheatgrass (Pseudoroegneria spicata Pursh A. Löve), Sandberg's bluegrass (Poa secunda J. Presl), and big sagebrush (Artemisia tridentata Nutt.) (Table 1).The Middle Grassland Zone has a higher cover of bluebunch wheatgrass and less bare ground than the Lower zone and big sagebrush is thought to represent a seral stage as a result of livestock grazing.Rough fescue (Festuca campestris Rydb) is the dominant plant species in the Upper Grassland Zone.
Soil development across Lac du Bois Grassland has been heavily influenced by climate and the mix of postglacial sediment-surface forms dominated by tills with hummocky surfaces overlaid by glacio-lacustrine silt and sand deposits at lower elevations (Watson 1977;van Ryswyk and McLean 1989).There has been a range in elevation boundaries between grassland zones reported at Lac du Bois (e.g., Spilsbury and Tisdale 1944;Tisdale 1947;van Ryswyk et al. 1966;McLean and Marchand 1968) due to uneven topography, variable aspect, and localized vegetation changes over time (Nicholson et al. 1991).In this study, boundaries were defined by Chernozemic great groups according to Munsell soil colour values of the Ah horizon using criteria outlined by the Canadian Soil Classification Working Group (SCWG 1998).Accordingly, Brown Chernozems were sampled from 389 to 628 m asl, Dark Brown Chernozems between 628 and 803 m asl, and Black Chernozems between 803 and 981 m asl (Fig. 1).These elevation bands for Brown, Dark Brown, and Black Chernozems align closely with the three distinct veg- ¶ Herbaceous aboveground biomass (excluding big sagebrush) was determined on a dry weight basis over 4 years (2017)(2018)(2019)(2020) for each grassland community by clipping green standing biomass as close to the soil surface as possible within four 1 m 2 frames at peak biomass (June-July) and drying at 60 • C for 48 h.
etation zones occurring at Lac du Bois, otherwise known as the Lower, Middle, and Upper Grassland Zones, respectively (Lee et al. 2014).

Microclimate across elevation gradient
Monthly and annual climate normal data (1991-2020) for all sampled positions (n = 24) along the elevation gradient (Fig. 1) were generated with ClimateBC v7.10 software package (2021) using downscaled PRISM climate data as described by Wang et al. (2016) and minimum and maximum values were grouped among soil zones (Table 2).At the lowest elevation of the gradient (389 m asl), mean annual precipitation (MAP) is 295 mm and mean annual temperature (MAT) is 9.0 • C. The MAP increases to 411 mm and MAT decreases by 3.5 • C at the highest point along the gradient (981 m asl).The sequence of soil and plant community zones across the elevation gradient is strongly controlled by soil water availability, which increases along the gradient and can be evaluated by the range in annual heat moisture (AHM) index that is greatest at lowest elevation (64.3) and decreases to 38.6 at the top of the gradient (Fig. 1).Potential evapotranspiration was calculated using the Hargreaves' equation (Hargreaves 1994) and was estimated at 734 and 654 mm for the lowest (389 m asl) and highest (981 m asl) points along the gradient, respectively.Soil development and associated plant communities are also affected by the number of frost-free days, which are estimated to be 234 days at the lowest elevation and 187 days at the highest elevation along the gradient over the last 30 years.A summary of climate data for the three soil zones is given in Table 2.

Selection of soil sample locations
Sampling site locations were selected within 25 m wide bands along the 600 m elevation gradient for a total of 24 sites from 400 to 1000 m asl (Supplemental S1).As a result, there were a total of nine sites in the Brown Chernozem zone, seven in the Dark Brown, and eight in the Black Chernozem zone (Fig. 1).Because of the complex topographical features of the study area, we followed the same approach for site selection as used in the 1961 study by van Ryswyk et al. (1966): (i) slope 10%-20%, (ii) slope length of 30 to 60 m, (iii) no curvature, and (iv) south-or southwest-facing aspect.In spring of 2019, potential sampling sites that met our topographic criteria were identified and mapped within each elevation band (n = 24) from a LiDAR derived digital elevation model with a 3 m × 3 m cell size.Potential sampling sites were then located in the field with a RTK GPS system (accuracy <1 m).Final site selection was based on local terrain features (e.g., plant species composition, avoidance of rock outcrops, old roadways, and obvious signs of human or livestock disturbance), as well as considerations of vehicle and hiking access.Our selection criteria were similar to those used when selecting "zonal" sites for BEC classification where the mature or climax plant community reflects the average climate for a particular biogeoclimatic subzones, rather than unique topography and soil conditions (MacKenzie and Meidinger 2018; Ryan et al. 2021).

Soil sampling and bulk density determination
Soil sampling (at the 0-5, 5-15, and 15-30 cm depths) occurred at each of the 24 elevation bands across the elevation gradient between 1 and 28 May 2019, using undisturbed soil cores (Blake and Hartge 1986).After careful removal of plant litter and cryptogamic crust, surface soil (0-5 cm) was sampled using a composite of five 35.5 cm 3 cores sampled with a drop hammer within a 30 × 30 cm area for a total sample volume of 177.5 cm 3 .Soil from 5-15 cm and 15-30 cm depths were sampled in the middle of each depth interval using one double-cylinder, drop hammer sampler, with a 7.5 cm diameter by 7.5 cm deep core (331.35cm 3 ).At each elevation band (n = 24), two profiles were sampled 5 to 10 m apart for a total of 48 unique soil profiles characterizing the 0-30 cm soil depth.Soil bulk density was calculated on a coarse fragmentfree basis as the mass of dry soil per volume of field-moist soil (Blake and Hartge 1986).The volume of coarse fragments was determined using a particle density of 2.65 g cm −3 and their air-dry mass.
In the spring of 2021, a soil auger (8.3 cm diameter) was used to sample the soil adjacent to the 2019 locations from each elevation band (n = 24) from 0-15, 15-30, 30-45, and 45-60 cm soil depths for further chemical and particle size analysis.Soil bulk density and coarse fragment estimates for the 30-45 and 45-60 cm soil depths were sampled from four pits in each soil zone using the excavation method as described by Page-Dumroese et al. (1999) and Throop et al. (2012).All soil  Note: AHM, annual heat moisture index; MAT, mean annual temperature; MAP, mean annual precipitation; PET, potential evapotranspiration calculated using Hargreaves's equation (Hargreaves 1994); FFD, frost-free days.
samples were air-dried to constant mass and passed through a 2 mm sieve to determine coarse fragments (>2 mm) before being shipped for analysis at the Analytical Chemistry Services Laboratory at the Ministry of Environment and Climate Change Strategy, Victoria, BC.

Soil chemical properties
Soil pH, total soil C, and SIC were determined on soil samples collected in 2019 and 2021.Prior to analysis, samples were ring ground to <1 mm and total soil C was determined by dry combustion (Nelson and Sommers 1996) using a Flash 2000 elemental analyser (Thermo Fisher Scientific, Waltham, MA, USA).Soil inorganic C was determined using a Primacs SNC-100 analyzer (Skalar, The Netherlands) that acidifies the sample using phosphoric acid and collected CO 2 is then measured using an IR detector.Soil organic C was calculated as the difference between total soil C and SIC.Soil pH was measured in a 0.01 mol L −1 CaCl 2 1:2 solution according to Hendershot et al. (2008).

Soil texture analysis
Soil texture was determined in 2021, on 96 unique samples taken from the 24 elevation bands over four depths (0-15, 15-30, 30-45, and 45-60 cm).Soil texture was determined by the integral suspension pressure (ISP+) method with a PARIO™ system (METER Group AG, Munich, Germany) as de-scribed by Durner et al. (2017) and Durner and Iden (2021).The ISP+ instrument and software uses an automated technique to determine size distribution of soil particle using an inverse modelling approach of suspension pressure data in a 1 L sedimentation cylinder.Before ISP+ analysis, soil samples are dispersed using the methodology as outlined for the hydrometer method (Gee and Bauder 1986).Soil samples (50 g) without any pretreatment for organic matter or carbonates were dispersed overnight using 100 mL Na 6 (PO 3 ) 6 (hexametaphosphate 50 g L −1 ) and 250 mL distilled H 2 O and then agitated on a mixer for 5 min.Dispersed samples were then transferred to sedimentation tubes and filled to 1000 mL with distilled H 2 O.

Carbon stock calculations
Carbon stocks (kg of SOC m −2 ) were determined for each sampled soil layer and summed to create the 0-30 and 0-60 cm stock estimates using eq.1: where SOC (g g −1 ) refers to the concentration of SOC, ffBD (kg m −3 ) the fine fraction (<2 mm) soil bulk density, and CF (m 3 m −3 ) the volume of coarse fragments (>2 mm).The thickness of the soil depths used in this study (0-5, 5-15, 15-30, 30-45, and 45-60 cm) reflects the standard sampling depths Table 3 Soil properties of the three Chernozemic great groups (here referred as zones) along the elevation gradient at Lac du Bois Grassland.as suggested by the FAO soil C sampling protocol (Yigini et al. 2018).To compare C stocks between 2019 and 1961 sampling events, we modified the 0-30 cm C stock calculation to match the 0-20 cm depth used by van Ryswyk et al. (1966) by using a depth of 5 cm, rather than 15 cm, from the 15-30 cm C stock calculation and then summed the three layers (0-5, 5-15, and 15-20 cm).To calculate the 0-60 cm soil C stock, we used the C concentration from the 0-30 cm depths sampled and analyzed in 2019 and then the 30-60 depth was calculated using the data obtained in 2021.

Statistical analysis
Mean values and standard deviation were used to characterize soil properties across Brown (n = 9), Dark Brown (n = 7), and Black (n = 8) soil zones along the elevation gradient.Linear and curvilinear regression analyses were performed in Sigmaplot 14.5 (Systat Software, San Jose, CA, USA) following a test for normality (Shapiro-Wilk) and constant variance (Spearman's rank correlation) among the residuals.Coefficient of determination (R 2 ) and standard error of the estimate (SEE) were used to select best fit between linear and curvilinear models.

Soil properties occurring along the elevation gradient
Across all sample locations, there was a similar trend of increasing sand and clay and a decrease in silt with depth (Table 3 and Supplemental S2).Similar to what was reported by van Ryswyk et al. (1966), there were only subtle changes in soil texture across the three soil zones, where silt-loam textures were most common in the Brown and Dark Brown zones and sandy-loam in the Black Chernozems.Within the Brown Chernozem zone, soil pH was 6.8 ± 0.3 at the surface (0-15 cm) and 8.0 ± 0.2 between 45 and 60 cm (Table 3).Soil pH in the Dark Brown Chernozem zone was 6.7 ± 0.2 and 8.0 ± 0.1 from the 0-15 and 45-60 cm soil depths, respectively.Black Chernozems had a pH of 6.5 ± 0.3 at the surface (0-15 cm) and 7.6 ± 0.5 from the 45-60 cm depth.Overall, there was a common trend of increasing pH with depth (6.7 ± 0.3 to 7.9 ± 0.3) and of decreasing pH with increasing elevation (7.5 ± 0.6 in the Brown zone, 7.4 ± 0.6 in the Dark Brown, and 7.1 ± 0.7 in the Black soil zone).
Soil bulk density at the soil surface (0-5 cm) decreased across the elevation gradient from 1245 ± 151 kg m −3 in the Brown zone to 1034 ± 100 kg m −3 in the Dark Brown zone and 771 ± 214 kg m −3 in the Black soil zone (Fig. 2B).A similar decrease across soil zones was measured at the 5-15 cm soil depth, with soil bulk density of 1161 ± 80 kg m −3 in the Brown zone, 983 ± 84 kg m −3 in the Dark Brown zone, and 878 ± 134 kg m −3 in the Black Chernozem.Soil bulk density from the 15-30 cm depth had the same negative relationship with increasing elevation, but the decrease went from 1100 ± 156 to 904 ± 130 kg m −3 along the Brown to the Black Chernozem sequence.Negative linear relationships among the three soil depths explained 58% (0-5 cm), 65% (5-15 cm), and 25% (15-30 cm) of the variability in soil bulk density.Soil bulk density from the 0-30 cm depth, calculated as a weighted average to reflect the different soil thickness of each depth measurement, showed a consistent negative relationship with increasing SOC across the Chernozem sequence with values of 1144 ± 147 kg m −3 (Brown), 1008 ± 153 kg m −3 (Dark Brown), and 873 ± 173 kg m −3 (Black).
Across all soil zones, coarse fragment (>2 mm) content (vol.%) at the soil surface (0-5 cm) was, on average, 2% and was never greater than 6% (Fig. 2C).Coarse fragments increased slightly at the 5-15 cm depth from a mean ± SD of 7 ± 12% in the Brown to 16 ± 16% in the Dark Brown and 13 ± 10% in Black Chernozems.Finally, the volume of coarse fragments was highest from the 15-30 cm soil depth with 14 ± 14% in the Brown soil zone, 28 ± 13% in the Dark Brown zone, and 22 ± 10% in the Black soil zone.A previous study by Evans et al. (2012) done in the Brown and Dark Brown soil zones at Lac du Bois reported similar coarse fragment volumes from comparable soil depths with 7 ± 7% (0-7.5 cm), 14 ± 9% (7.5-15 cm), and 18 ± 8% (15-30 cm).Increasing coarse fragment content with depth, while also increasing with elevation, reflects the postglacial deposit of an aeolian veneer up to 10 cm thick that blanketed the morainal sediments with fine sand and silt deposits following the rapid drainage of Glacial Lake Thompson (Watson 1977;van Ryswyk and McLean 1989;Johnsen and Brennand 2004).Glacial Lake Thompson has been estimated to cover much of the Brown and some of the Dark Brown soil zones before draining, leaving large deposits of silt and fine sand that are still prominent today in much of the Thompson River Valley.
SOC stocks (kg m −2 ; 0-30 cm) increased from an average of 3.64 ± 1.0 in the Brown zone to 3.95 ± 0.73 in the Dark Brown and 6.77 ± 1.4 kg m −2 in the Black Chernozem zones (Fig. 2D).A curvilinear model best fit the data across the elevation gradient (Fig. 2D), explaining 83% of the variability compared with 54% explained by the linear model.The nonlinear relationship was driven by the high variability in SOC (g kg −1 ; Fig. 2A) within the Brown soil zone.A long-term grazing trial located within the Dark Brown soil zones at the Lac du Bois Grassland (Evans et al. 2012) reported total C stocks (0-30 cm) with a range of 5.00 to 5.99 kg m −2 , which are higher than our estimates for Dark Brown Chernozems, but are likely a reflection of the range in SOC from different landscape positions within this zone.
C stock comparison with data reported by van Ryswyk et al. (1966) Original soil C stocks data from the elevation gradient (van Ryswyk et al. 1966) were extracted using WebPlotDigitizer v4.5 software and means were summarized for the 0-20 cm soil depth for Brown (2.53 ± 0.43 kg m −2 ), Dark Brown (4.17 ± 0.73 kg m −2 ), and Black (7.89 ± 1.8 kg m −2 ) soil zones.After 58 years (Fig. 3A), SOC stocks (kg m −2 ) in the Brown soil zone were similar to the original measurements (2.54 ± 0.62) but were reduced by 26% across Dark Brown Chernozems (4.17 to 3.1 kg m −2 ), and by nearly 35% for Black Chernozems (7.89 to 5.14 kg m −2 ).The parameters used in the original SOC calculation (van Ryswyk et al. 1966) were vague with respect to how mineral particles >2 mm were accounted for.Hence, we calculated SOC stocks (Fig. 3) using three different equations as described by Poeplau et al. (2017).In Fig. 3A, we used fine fraction bulk density (ffBD) and an additional parameter (1 − CF m 3 /m 3 ), described as Method 4 (M4) in Poeplau et al. (2017), so that mineral particles <2 mm were accurately accounted for in the SOC stock estimates.In Fig. 3B, we used ffBD only (Poeplau et al.'s Method 2,M2), producing SOC estimates that were 10% greater in Brown Chernozems, 13% smaller in Dark Brown Chernozems, and 27% smaller in Black Chernozems when compared with the 1961 measurements.In Fig. 3C, no corrections for coarse fragments were made in either the mass or volume of the sample (Poeplau et al.'s Method 1, M1), producing SOC estimates that were 20% Fig. 3. Soil organic carbon (SOC; kg m −2 ; 0-20 cm) across the elevation gradient in 2019 and 1961, where the latter were taken from the study by van Ryswyk et al. (1966).Using a factor of 10, we converted the original SOC unit (g 2 dm −3 ) to kg m −2 (1 m 2 = 100 dm 2 and 1 kg = 1000 g).Note that a similar sampling approach and design were used in 1961 and 2019, but three different SOC stock calculations are used in the 2019 data as described by Poeplau et al. (2017): (A) M4 calculates soil C stocks with coarse fragments in the soil volume and the fine fraction soil mass; (B) M2 calculates soil C stocks with no coarse fragments (>2 mm) in the soil volume and the fine fraction soil mass; and (C) M1 calculates soil C stocks using a bulk density value not corrected for coarse fragments.and 7% higher in the Brown and Dark Brown zones, respectively, but 12% lower in the Black Chernozems compared with the 1961 measurements.This study highlights why SOC stock estimates from rocky soils are challenging (Page-Dumroese et al. 1999;Throop et al. 2012;Smith et al. 2020) and selecting the wrong equation can significantly over-or underestimate SOC stocks, thereby confounding future comparisons (Poeplau et al. 2017).
Selecting an appropriate method to calculate SOC stocks is not the only challenge when evaluating the effects of grazing regimes on soil C sequestration.Livestock grazing can also affect SOC contents due to a change in soil compaction (e.g., Naeth et al. 1990;Henderson et al. 2004;Abdalla et al. 2018).Past measurements of soil penetration resistance and bulk density in the Black Chernozem zone (Krzic et al. 2014) and the Brown and Dark Brown soil zones (Evans et al. 2012) at Lac du Bois Grassland have a decrease in these physsoil properties with grazing exclusion (20 to 60 years).Those findings demonstrate the impact that livestock grazing can have on soil compaction and subsequent estimates of SOC stocks using soil bulk density.When comparing C stocks among different levels of grazing intensities, it has been suggested that a C mass per equivalent soil mass approach should be used (Henderson et al. 2004) to combat the confounding effect that grazing animals have on soil bulk density and the subsequent SOC stock calculations (Ellert and Bettany 1995;von Haden et al. 2020).Sampling over multiple depths and providing the proportion of fine (<2 mm) and coarse fragments (>2 mm) within the sampling depth are important for SOC estimates.This is of particular importance for studies that elect to use an equivalent soil mass, rather than a fixed soil depth approach, to evaluate SOC stocks change with grazing management.
Similar to grazing-induced changes to soil bulk density, grazing can also modify plant community composition.The reduction in Kentucky bluegrass (Poa pratensis L.) cover in the Upper Grassland from 42% in 1961 (van Ryswyk et al. 1966) to 12% in this study (Table 1) could help to explain the surprisingly high SOC stocks from the 0-20 cm soil depth in the original study.The introduction and prevalence of Kentucky bluegrass, an exotic species with higher root density at the soil surface, has been associated with increased SOC (Sanderson et al. 2017;Bork et al. 2020) in part due to a buildup of dead biomass (thatch) at the soil surface with lower C:N ratio (Toledo et al. 2014).Spreading by rhizomes, Kentucky bluegrass will also form a continuous mulch on the surface, which has been shown to regulate soil water and temperature dynamics (Avery et al. 2019), reduce plant diversity and germination (Halvorson et al. 2022), and lead to increased soil C and N (Sanderson et al. 2017;Hendrickson et al. 2021).Despite the documented decrease in Kentucky bluegrass from the Upper Grassland and Black Chernozem zone, this species has been increasing on mesic grasslands with and without grazing in both Alberta (Zapisocki et al. 2022) and BC (Bradfield et al. 2021;Hamilton et al. 2022), and throughout the Northern Great Plains (Toledo et al. 2014).Where present, the influence of Kentucky bluegrass needs to be taken into consideration in future studies on grassland SOC change and storage.
Because we followed the relatively strict topographical features for site selection outlined by van Ryswyk et al. (1966), we are satisfied that our resampling on the elevation gradient at Lac du Bois Grassland after 58 years was an accurate comparison, despite not knowing the exact locations used in the original study.In both studies, SOC increased more rapidly with elevation at or near the 600 m asl in the Dark Brown and Black soil zones.The influence of reduced water deficits or decreasing AHM (Fig. 1) with increasing elevation on SOC was similar to that first described van Ryswyk et al. (1966).Numerous studies have shown that climate and other abiotic factors have a stronger influence on SOC dynamics than management alone (Brown et al. 2010;Augustin and Cihacek 2016;Deane McKenna et al. 2022), which may be a reasonable explanation for the lack of SOC change with increasing C inputs.Our improved understanding on the influence of elevation and microclimate on the distribution of SOC along this elevation gradient will allow for future quantifications of management effects on soil C sequestration.
Total, organic, and inorganic C stocks within the 0-60 cm soil depth Bulk density measurements used in C stock calculations from the 30-45 cm soil depth (data not shown) were 1400 kg m −3 in Brown, 1180 kg m −3 in Dark Brown, and 1200 kg m −3 in Black soil zones.Soil bulk density from the 45-60 cm depth was 1227 kg m −3 in the Brown soil zone, 1049 kg m −3 in the Dark Brown zone, and 1175 kg m −3 in the Black Chernozem.The volume of coarse fragments from the 30-45 cm depth was 37% in the Brown and Dark Brown zones and 30% in the Black soil zone.Similarly, coarse fragments from the 45-60 cm soil depth were 43% in the Brown soil zone and 30% in both the Dark Brown and Black soil zones.
Total C stocks (Fig. 4A; 0-60 cm) increased from 7.48 ± 2.1 kg m −2 in the Brown soil zone to 8.59 ± 1.4 kg m −2 in the Dark Brown zone, and 12.37 ± 3.8 kg m −2 in the Black soil zone.A positive linear relationship explained 45% of the variation in total soil C with elevation compared with the curvilinear model (R 2 = 0.69).SOC stocks (kg m −2 ) to a depth of 60 cm (Fig. 4B) were 5.73 ± 1.7 in Brown, 5.87 ± 0.76 in Dark Brown, and 10.11 ± 2.5 in Black Chernozems.In comparison, estimates of SOC stocks from Canadian Prairie soils ranged from 8.0 to 15.0 kg C m −2 across Brown, Dark Brown, and Black Chernozems (Landi et al. 2003a;Pennock et al. 2011).Those estimates of SOC stocks from Canadian Prairies were higher since they were obtained from deeper soils with less coarse fragments than our estimates at Lac du Bois.
Soil inorganic C stocks (Fig. 4C) represented an average of 23% of the total C stock in the Brown zone and 32% and 18% from the Dark Brown and Black soil zones, respectively.Depth to SIC also increased with elevation where 80% of SIC was measured below the 30 cm depth in the Brown soil zone and 88% and 99% in Dark Brown and Black Chernozems, re-spectively.Total C stocks (Mg ha −1 ; 0-60 cm) in Brown and Dark Brown Chernozems from Saskatchewan (Bork et al. 2020) were nearly twice those measured in our study at Lac du Bois, yet SIC measured in their study represented a slightly higher proportion (19% to 47%) of the total soil C mass.Omitting SIC from SOC calculations may lead to an overestimation of SOC stocks, especially if deeper soil C stocks are to be considered on grassland soils that are high in carbonates (Reeder and Schuman 2002).
The greater SOC variability in the Brown soil zone was due to the largest C stock being measured at the lowest elevation in the Brown Chernozem zone (circled 391 m asl; Fig. 4).This location had a SOC stock (10.01kg m −2 ) that was nearly twice the average (5.2 kg m −2 ) of all other elevations within the Brown soil zone.The linear model of SOC with elevation was greatly improved (R 2 = 0.69) when this data point was removed from the analysis (data not shown) with the corresponding equation of y = 0.012x − 1.666.Likewise, the polynomial equation explained 84% of the SOC variability with elevation (y = −0.044x+ 3.95e−5x 2 + 17.08) when this point was excluded.The high soil C stock (0-60 cm) from the lowest elevation along the gradient was not caused by an extremely high concentration of C at any particular depth, but due to the uniform distribution of soil C within each of the four soil depths sampled (Supplemental S3).Soil texture measured from these depths (Supplemental S2) also indicated a uniform distribution with depth with virtually no clay measured from the 0-60 cm profile and an almost even proportion of sand and silt across the four depths that ranged from 50%-63% for sand and 37%-50% for silt.Considering the consistent soil C concentration and relatively uniform texture at all four sampling depths, it is likely that soil at this sampling location reflects buried Ah horizons on alluvial fans formed from outwash material following one or several destabilization events since glaciation (e.g., Pennock et al. 2011;Helgason et al. 2014).Recognizing these alluvial fans and outwash material is important when selecting sample sites in grassland ecosystems with complex and steep terrain.5-15, and 15-30 cm) and 2021 (0-15, 15-30, 30-45, and 45-60 cm).(B) Soil pH vs. ratio of SIC to total C at Lac du Bois Grassland, north of Kamloops, British Columbia.
Soil inorganic carbon and soil pH SIC had a noticeable relationship with soil pH, where concentrations of SIC above detectable limits (0.03%) were most common when soil pH > 7.0 (Fig. 5B).Soil samples that had >50% of total C composed of SIC were exclusively sampled from the 30-60 cm soil depths and had soil pH values higher than 7.5.Although van Ryswyk et al. (1966) did describe carbonate-enriched C horizons (Ck) in all Chernozem profiles, no estimates of SIC concentrations with depth were given in their study.
In the conterminous US, as much as 40% of all soil C is composed of SIC, but land-use effects on SIC stocks are not well understood compared with the effects on the SOC cycle (Lorenz and Lal 2022).Under alkaline conditions, dissolved CO 2 from increased soil respiration or elevated atmospheric concentrations and an abundance Ca 2+ in soil solution can promote the precipitation of CaCO 3 and serve as a pathway for atmospheric CO 2 sequestration into the SIC pool (Reeder et al. 2004;Zamanian et al. 2016).When compared with the SOC cycle, SIC is thought to be somewhat insensitive to the effects of climate change, such as rising temperatures and increased variability of incoming precipitation, and therefore has much longer residence times in soil (Heckman and Rasmussen 2018).Furthermore, the rate of atmospheric CO 2 removal has been estimated to be greater from the SIC cycle than what can be sequestered by the SOC cycle (Landi et al. 2003b).Hence, there is equal value in recognizing the C sequestration potential of the SIC cycle in studies from semiarid grassland ecosystems (Reeder and Schuman 2002;Zamanian et al. 2016), rather than only focusing on SOC.

Conclusions
Consistent with the previous study from 1961 carried out on this elevational gradient (van Ryswyk et al. 1966), SOC stocks (0-60 cm) increased from 5.73 kg m −2 in the Brown zone to 5.87 kg m −2 in the Dark Brown and 10.11 kg m −2 in the Black Chernozem zones.Soil inorganic C represented 23%, 32%, and 18% of the total C stock among Brown, Dark Brown, and Black soil zones, respectively, with the highest concentrations located below 30 cm and soil pH > 7. Despite the improved grazing regime since the original soil C measurements in 1961, there was not a consistent increase or decrease in the 0 to 20 cm soil C stock across the three soil zones.The lack of clear information on the parameters used to calculate SOC stocks by van Ryswyk et al. (1966) and the known biases that occur with discrepancies in methodologies used to calculate SOC stocks in rocky soils make the comparison over time problematic.In the Brown soil zone, SOC stocks increased by 10%-20%, whereas in the Black soil zone SOC decreased by as much as 35%.Reductions in Kentucky bluegrass and an increase in deeper rooted bunchgrass species in the Black Chernozem zone may have had a greater effect on SOC stocks below the 0-20 cm depth, which could explain the decrease in SOC in this zone.The effect of improved grazing management on soil C sequestration over the last 58 years remains unclear at Lac du Bois Grassland.However, the positive exponential relationship that SOC has with elevation, as first described by van Ryswyk et al. (1966), was evident along the Chernozem sequence from all soil depths, which has provided a better understanding of how elevation and climate influence the distribution of Chernozemic soils across this landscape.Future assessments of management-and climateinduced changes on SOC along the elevation gradient will benefit from our new measurements that document both SIC and coarse fragments at greater depths.

Fig. 1 .
Fig. 1. (A) Elevation boundaries for the Brown, Dark Brown, and Black Chernozem soil zones and associated plant communities at the Lac du Bois Grassland, north of Kamloops, British Columbia.(B) Elevation and calculated annual heat moisture index (AHM) of sample locations used in the evaluation of soil carbon stocks among Chernozemic great groups (Brown, n = 9; Dark Brown, n = 7; Black, n = 8).

Fig. 4 .
Fig. 4. Total C stock (kg m −2 ) (A), organic C stock (kg m −2 ) (B), and inorganic C stock (kg m −2 ) (C) to 60 cm depth across the elevation gradient at Lac du Bois Grassland, north of Kamloops, British Columbia.Alternate regression analyses were completed following the removal of the circled point at 400 m above sea level and discussed in text.

Table 1 .
Ground and vegetation cover (%) from the Lower, Middle, and Upper Grassland Zones and aboveground biomass at Lac du Bois Grassland.

Table 2 .
Modelled climate variables (ClimateBC) of the three Chernozemic great groups (here referred as zones) along the elevation gradient at Lac du Bois Grassland, north of Kamloops, British Columbia.