Vertical Distribution of Soil Organic Carbon Density in Relation to Land Use / Cover , Altitude and Slope Aspect in the Eastern Himalayas

In-depth understanding about the vertical distribution of soil organic carbon (SOC) density is crucial for carbon (C) accounting, C budgeting and designing appropriate C sequestration strategies. We examined the vertical distribution of SOC density under different land use/land cover (LULC) types, altitudinal zones and aspect directions in a montane ecosystem of Bhutan. Sampling sites were located using conditioned Latin hypercube sampling (cLHS) scheme. Soils were sampled based on genetic horizons. An equal-area spline function was fitted to interpolate the target values to predetermined depths. Linear mixed model was fitted followed by mean separation tests. The results show some significant effects of LULC, altitudinal zone and slope aspect on the vertical distribution of SOC density in the profiles. Based on the proportion of mean SOC density in the first 20 cm relative to the cumulative mean SOC density in the top meter, the SOC density under agricultural lands (34%) was more homogeneously distributed down the profiles than forests (39%), grasslands (59%) and shrublands (43%). Similarly, the SOC density under 3500–4000 m zone (35%) was more uniformly distributed compared to 3000–3500 m zone (43%) and 1769–2500 m and 2500–3000 m zones (41% each). Under different aspect directions, the north and east-facing slopes (38% each) had more uniform distribution of SOC density than south (40%) and west-facing slopes (49%). OPEN ACCESS


Introduction
Soil, as the largest terrestrial carbon (C) sink, stores almost three times the amount of C in the vegetation and approximately double that in the atmosphere [1].Globally, the soil stores between 1500 and 1600 Pg•C in the upper meter [2][3][4] and about 2350 Pg•C in the upper 3 m [5].The amount of stored soil organic carbon (SOC) is the net balance between C inputs from plant production and losses through decomposition and leaching [6].Hence, depending upon C fluxes, the soil can either be a source or sink for atmospheric CO 2 [7,8].The stability and distribution of SOC in the soil profile is influenced by biotic controls such as the abundance and vigor of faunal, microbial and plant species as well as environmental controls like temperature, moisture and soil texture [9].
Jobbàgy and Jackson [5] reported a strong relationship between vertical distribution of SOC and climate, vegetation and soil texture with an estimated 1502, 491 and 351 Pg•C for the first, second and third meter, respectively.Climate exerts a first order control on the amount [3] and cycling rate [10] of SOC at the global scale.Temperature, moisture (precipitation) and solar radiation are the most important factors influencing plant growth, litter production and SOC mineralization [11].While SOC content increases with precipitation, it decreases with temperature due to higher decomposition of organic matter (OM).However, the sensitivity of SOC decomposition to temperature becomes decreasingly acute at higher than at lower temperatures [12].This is critical for maintaining high SOC stocks in high-altitude areas, such as the Himalayas, in response to global warming.
Topography (especially altitude and slope aspect) through its influence on climate, soil properties and soil water content [11], is largely responsible for the spatial and vertical distributions of SOC.The expected decrease in temperature with altitude depresses OM decomposition rates more than litter production, and generally favors the accumulation of SOC.However, as the soil becomes shallower at the higher altitudinal zones, the overall SOC storage capacity is limited despite higher SOC content per unit of soil mass.While a number of previous studies reported increased SOC with altitude [13,14], there have been records of decreases in SOC with altitude [15,16].
Because of these complex environmental interrelationships, SOC is highly heterogeneous in mountainous areas mainly as a result of local-scale variability in soil environment and microclimate [17,18].The strong relationship of slope aspect with SOC is through the former controlling the microclimatic factors such as soil temperature, moisture, vegetation and microorganisms by affecting the solar radiation and evapotranspiration.While many studies have reported higher SOC content [19][20][21] on the north-facing slopes of the Himalayas, Sidari et al. [18] reported lower SOC content on the northern aspect direction.In contrast, Han et al. [22] did not find any significant difference in SOC content between the north and south-facing slopes.
Similar to altitude and aspect, vegetation cover and land use have significant impact on SOC dynamics through organic carbon (OC) inputs, decomposition and stabilization.The plant functional type, which is specific to a particular plant species, has a significant effect on the vertical distribution of SOC [5].The vegetation cover not only determines the amount and type of OC inputs but also their Land 2014, 3 1234 location, as either aboveground litter deposited on the soil surface or in the subsoil as root biomass.It also influences the decomposition and stabilization of SOC in the soil [9].High-altitude ecosystems are reported to have generally high root: shoot ratios than other ecosystems [23,24].This leads to an expectation that SOC distribution with depth is shallower in high-altitude areas.While this might provide a unique vertical distribution of SOC, very little is currently known about it [5,9,25].
Thus, characterization of vertical distribution of SOC density is important for C budgeting, designing appropriate C sequestration strategies and enhancing our understanding of biogeochemical cycles in different ecosystems.It is also crucial for examining the responses of terrestrial ecosystems to global change through biogeochemical modelling [5,25,26].However, such studies are very few at the regional [25], national [26] and global [5] scales, especially for high-altitude ecosystems such as the Himalayas.The aim of this study was therefore to investigate how the vertical distribution patterns of SOC density in the profiles varied under different land use/land cover (LULC) types, altitudinal zones and aspect directions in a montane ecosystem of Bhutan, Eastern Himalayas.

Study Area
The study was done in the northern part of Wang Watershed in Western Bhutan which is located between 89°14.03′E to 89°35.82′E and 27°04.87′N to 27°38.28′N (Figure 1).The study area (1014 km 2 ) as described in Dorji et al. [27] represents a montane ecosystem in the Himalayas with high peaks and deep valleys.The altitude ranges from 1769 to 5520 m above sea level (a.s.l) with climatic conditions varying from warm temperate in the valley floor to cold and dry alpine on the ridges.
While the mean maximum annual temperature for the 1769-2500 m zone is about 20 °C, the mean minimum annual temperature is 8 °C (Table 1).Similarly, the 2500-3000 m zone has a mean maximum annual temperature of approximately 17 °C and the mean minimum annual temperature of about 7 °C (Table 1).The zone above 3000 m expectedly experiences lower mean annual temperature as the lapse rate is about 5 °C per vertical 1000 m.Rainfall is monsoonal with mean annual precipitation varying from about 755 mm in low altitudinal zone (1769-2500 m) to about 2990 mm in the mid-altitudinal zone (2500-3000 m) (Table 1).Although there is no rainfall data available for the zone above 3000 m, the thick forest cover between 3000 and 4000 m is indicative of higher mean annual rainfall up to 4000 m, which probably slowly declines in the upper altitudinal zone (4000-5520 m zone).Geologically, the study area is largely underlain by mixed grade metamorphic rocks [28,29].Detailed soil information of the study area is not available however weathered and leached thin dark topsoil over bright subsoil up to 3000 m, non-volcanic andosolic soils and weak podzols up to 4000 m and deep dark and friable topsoil over yellowish subsoil mixed with raw glacial deposits above 4000 m are dominant in Bhutan [30] (Table 1).Forests (73%) dominate the LULC of the study area followed by shrublands (14%), agricultural lands (7%) and grasslands (3%) [31].The detailed distribution of LULC types under different altitudinal zones is presented in Table 1.The study area was broadly grouped into five altitudinal zones (i.e., 1769-2500 m, 2500-3000 m, 3000-3500 m, 3500-4000 m and 4000-5520 m) and four aspect directions (i.e., north = 316°-45°, east = 46°-135°, south = 136°-225° and west = 226°-315°) to investigate how the vertical distribution of SOC density in the profiles varied under different altitudinal zones and aspect directions.

Sampling Strategy and Laboratory Analyses
The condition Latin hypercube sampling (cLHS) scheme [32] was used to locate the sampling sites using 90 m digital elevation model (DEM), normalized difference vegetation index (NDVI), SAGA wetness index (SWI) and LULC as covariates.Detailed descriptions about the derivation and extraction of these environmental covariates can be found in Dorji et al. [27].The cLHS scheme is based on a stratified random procedure, which provides an efficient way of sampling from the multivariate distribution of the covariates.It identifies a set of values from several covariates that satisfies the requirement of Latin hypercube of only one sample in each row and column in n dimensions.Since it is cumbersome to select the best combinations of values that exist in the covariate data, many iterations are required to establish suitable combinations.However, the iterations can be minimized by selecting samples that form a Latin hypercube in the feature space [32] such that from a given sites with ancillary data ( ), it selects sample sites ( ≪ ) and the sampled sites form a Latin hypercube [32].Thus, using the cLHS algorithm a total of 100 field sampling sites was generated covering all LULC types, altitudinal zones and aspect directions (Figure 1).Of the 100 sampling points, only 93 sites could be sampled due to poor accessibility.Paired sampling was done to account for the local variability of the target variable by locating each paired site at a random distance of 1-2 km from the primary site.At each sampling point, a soil profile pit was dug, described and sampled by genetic horizons to a depth of 1 m or bedrock.The presence of carbonates in the soil was tested using 1 M hydrochloric acid (HCl) in the field.Disturbed bulk density was measured by the core ring method [33].Buried organic materials (e.g., fine roots and dead wood fragments) were removed from the soil sample during air drying.The air-dried samples were sieved to 2 mm to carry out the chemical analyses.Soil pH was measured by pH meter in 1:2.5 suspension of soil in distilled water.The <2 mm air-dry soil samples were ground to <53 µm to determine the C content using a vario Max CNS Element Analyzer.The C content (OC) was corrected to an oven-dry basis.

Fitting Equal-Area Spline Profile Function
The horizon-based SOC concentration, bulk density and pH values were interpolated to predetermined depths of 0-20 cm, 20-40 cm, 40-60 cm, 60-80 cm and 80-100 cm by fitting an equal-area spline profile function using the CSIRO SplineTool V2 [34].The equal-area spline profile function is based on the assumption that all true soil attributes consistently vary with depth [35].This assumption can be represented mathematically by denoting depth by , and the depth function by ( ) such that, ( ) and its first derivative ′( ) are both continuous, and that ′( ) is square integrable.The depths of the boundaries of the layers are given by < , … < .Apart from the measurement error, the measurement of the bulk sample from horizon is assumed to reflect the mean attribute level.The measurements ( = 1, … ) are therefore mathematically modelled as where ̅ = ( ) /( ) is the mean value of ( ) over the interval ( , ).The errors are assumed independent, with mean 0 and common variance σ². ( ) represents a spline function and is determined by minimizing: The first and the second terms represent the fit to data and the roughness of function ( ) (expressed by ′( ) ), respectively.The trade-off between the fit and the roughness penalty is controlled by the parameter λ [35].

SOC Density Computation
The SOC density is C mass per unit area for a given soil depth and it was computed using SOC concentration, bulk density and depth: where = SOC density (kg•m −2 ), _ = SOC (kg/kg) in oven dry basis, = bulk density (kg•m −3 ) and = depth interval thickness (m) [27].In this study, coarse fragments (i.e., stone and gravels) were not accounted while computing the SOC density and this might potentially raise the SOC density values slightly.However, given that the coarse fragment content in the soil was generally low, exclusion is not expected to significantly alter results.

Data Analyses
A linear mixed model was used to analyze this data using SPSS Version 22.0 [36].The model was fitted to SOC density data at each depth interval using LULC type, altitudinal zone and slope aspect direction as fixed effects and sampling sites as random effect.Similarly, the model was also fitted for individual LULC type, altitudinal zone and aspect direction.A posteriori mean separation tests were performed using Fisher's LSD with Bonferroni adjustment [36].

Soil Properties under Different LULC Types
Table 2 shows the mean soil pH and bulk density values at different depths under different LULC types.The low pH values (<7.0) and the negative HCl test in the field confirmed absence or insignificant amount of carbonates in the soil and therefore all measured C was considered as OC.
The pH values at all depths were relatively low under forests (pH 5.3-5.8)compared to grasslands (pH 5.9-6.7),shrublands (pH 6.0-6.9) and agricultural lands (pH 6.2-6.9).This suggests that soils under forests are more acidic than other LULC types and this could be attributed to more leaching as a result of high precipitation under forests.On the other hand, the mean soil pH increased with depth under all LULC types.The mean bulk density values at all depths were higher under agricultural lands (1.11-1.34g•cm −3 ) than at the corresponding depths under forests (0.90-1.22 g•cm −3 ), grasslands (0.89-1.35 g•cm −3 ) and shrublands (0.99-1.23 g•cm −3 ).This is possibly due to low OM inputs and moderate soil compaction under agricultural lands due to tillage (ploughing and puddling) and livestock trampling.On the other hand, the relatively low bulk density under forests could be ascribed to high SOC and presence of non-volcanic andosolic soils characterized by moderately high SOC and very low bulk density [30,37].As expected, the mean bulk density under all LULC types slightly increased with soil depth (Table 2).

Influence of LULC, Altitudinal Zone and Aspect Direction on SOC Density
By combining SOC concentration, bulk density and soil depth thickness, we estimated the SOC density (Equation ( 3)).The results of the linear mixed model show some significant effects (p < 0.05) by LULC type, altitudinal zone and aspect direction on SOC density (Tables 3 and 4).In considering each soil depth separately, only altitudinal zone in the first three depths and aspect direction in the lowest depth had significant effects (p < 0.05) on SOC density (Table 3).Furthermore, the interactions between the fixed factors also had very limited impacts on SOC density except the influence of altitudinal zone × aspect direction interaction at 40-60 cm and 60-80 cm depths.This indicates that altitudinal zone has a stronger effect on SOC density compared to aspect direction and LULC.This is probably because SOC content generally increases with altitude due to high OM inputs and slow decomposition.The significant effect (p < 0.05) of the interaction between altitudinal zone × aspect direction on SOC density suggests that their combined effects influence SOC density in the deeper layers of the soil profiles.The effects of the fixed factors and their interactions on SOC density were much more significant when individual LULC type, altitudinal zone and aspect direction was considered (Table 4).The SOC density under all LULC types was significantly affected (p < 0.05) by both altitudinal zone and soil depth except under grasslands.Furthermore, the interactions between altitudinal zone × aspect direction under forests and shrublands, and altitudinal zone × depth under forests also had significant effects on SOC density.Similarly, the SOC density under all aspect direction was significantly influenced (p < 0.05) by LULC, altitudinal zone and depth except for north-facing slopes by LULC (Table 4).The interactions between LULC × altitudinal zone and altitudinal zone × depth also significantly affected (p < 0.05) the SOC density under east and south-facing slopes, and south-facing slopes, respectively.The SOC density under different altitudinal zones was however significantly influenced (p < 0.05) by aspect direction and soil depth except for the 4000-5520 m zone.Since none of the factors had significant effect on SOC density under 4000-5520 m zone, the results are not presented in Table 4.

SOC Density and Its Vertical Distribution under Different LULC Types
The mean SOC density values for the depths under different LULC types ranged from 2.5 to 6.0 kg•m −2 under agricultural lands, 2.5-9.6 kg•m −2 under forests, 1.1-11.3kg•m −2 under grasslands and 1.5-8.7 kg•m −2 under shrublands (Table 5).Except for the first depth interval, forests had relatively large mean SOC density than other LULC types.Mean SOC density was lowest under agricultural lands in the first three depths but was slightly larger than grasslands and shrublands in the last two depths.Although the differences in mean SOC density values among LULC types were clear, only SOC density under agricultural lands was significantly different (p < 0.05) from forests and grasslands at 0-20 cm depth and forests at 20-40 cm depth.Conversely, the cumulative mean SOC density in the top meter under different LULC types showed a similar trend with significantly larger (p < 0.05) cumulative mean SOC density under forests (24.84 kg•m −2 ) compared to shrublands (20.56 kg•m −2 ), grasslands (17.86 kg•m −2 ) and agricultural lands (17.38 kg•m −2 ).The cumulative mean SOC density in the first meter under forests was 7.5 kg•m −2 (43%), 7.0 kg•m −2 (39%) and 4.3 kg•m −2 (21%) larger than agricultural lands, grasslands and shrublands, respectively.The relatively large SOC density under forests and shrublands could be due to high OC inputs, slow decomposition and deep rooting depth compared to grasslands and agricultural lands.Furthermore, the presence of non-volcanic andosolic soils under forests, with moderately high SOC content throughout the solum [30,37]  While the mean SOC density invariably decreased with depth, its vertical distribution pattern was notably different under different LULC types.The mean SOC density under different LULC types ranged from 6.0-11.3kg•m −2 (0-20 cm), 3.9-6.4kg•m −2 (20-40 cm), 2.5-3.7 kg•m −2 (40-60 cm), 1.5-2.5 kg•m −2 (60-80 cm) and 1.1-2.5 kg•m −2 (80-100 cm) (Table 5).The mean SOC density values at 0-20 cm and 20-40 cm depths under all LULC types were significantly different (p < 0.05) from each other and also from the values of the last three depths.The exception is that the SOC density at 20-40 cm depth was not significantly different from the SOC density values at the last three depths under grasslands.On the other hand, the last three depths under all LULC types were not significantly different from each other, with the exception of values of 40-60 cm compared to those of the last two depths under forests.The relatively large mean SOC density in the top layers could be attributed to higher OM inputs by litter fall causing higher accumulation of SOC in the upper layer of the mineral soil.The effect of LULC on SOC density diminishes with depth and the vertical distribution pattern of SOC density largely depends on how deeply the OM is disseminated throughout the profile.
When considered in proportional terms, the percentage of mean SOC density in the top 20 cm relative to the cumulative mean SOC density in the first meter was highest under grasslands (59%) followed by shrublands (43%), forests (38%) and agricultural lands (34%) (Figure 2).These proportions indicate that more than 34% of the total SOC density in the upper meter is stored in the first 20 cm depth under all LULC types.The proportion of SOC density in the first 20 cm (relative to the top meter) also provides information about the vertical distribution of SOC density in the profiles.The higher the proportions of SOC density in the top 20 cm, the less homogeneous is the SOC density distributed down the profile.In this study, grasslands had the highest proportion of SOC density in the top 20 cm depth suggesting that SOC density is less homogeneously distributed down the profile compared to other LULC types.This is not unexpected as it could be attributed to very high OM inputs in the top layer which then decreases considerably down the profile.This could also be due to the very shallow rooting depth under grasslands [5,26].Comparatively, the proportion of SOC density in the first 20 cm (relative to the upper meter) was lowest under agricultural lands.This could be due to the opposite effect to that under grasslands and deep incorporation of farmyard manure and crop residue under agricultural lands.Although tillage is limited under orchards, it is moderately intensive under dry land agriculture and paddy land.Power tillers are often used at the study site for ploughing and puddling of paddy soils.Farmyard manure and/or leaf litters are applied every cropping season to improve the soil fertility.Conversely, the vertical distribution of SOC density under forests was more homogeneous than grassland and shrublands (Figure 2) probably due to high OM inputs from both above-and belowground biomass, slow decomposition and SOC complexation with andic and spodic properties of non-volcanic andosolic soils [37].Results reported here were, however, dissimilar to those from Wang et al. [26] in China and Jobbágy and Jackson [5] for the global ecosystems.The work by Wang et al. [26] reported highest proportion of SOC density in the first 20 cm (relative to the first meter) under forests (54%) followed by shrublands (46%), grasslands (39%) and cropland (37%).Similarly, Jobbàgy and Jackson [5] reported highest percentage under forests (50%) followed by grassland (42%), cropland (41%) and shrublands (33%).The abrupt decrease in SOC density with depth (relative to the SOC density in the first 20 cm) under grasslands (59%) reported here compared to the grassland (42%) reported for the global ecosystems [5] could be attributed to the shallower vertical root distribution under grasslands in the high-altitude ecosystems such as the Himalayas compared to the tropics.In support of this premise, a recent study by Yang et al. [24] found that 90% of the roots in the alpine grasslands were concentrated in the top 30 cm compared with 65% reported for the global ecosystems [23].Additionally, a recent study by Ota et al. [39] supports the influence of rooting-depth on the distribution of SOC in the soil profile.The SOC density under shrublands (43%) reported here was less homogeneously distributed down the profile compared to the shrublands (33%) reported for the global ecosystems [5] and this could also be due to shallower rooting depth under shrublands in the temperate region than in the tropics.However, the results for the shrublands reported here are comparable with the findings of a similar study reported by Wang et al. [26] in China.Under forests however, the SOC density (39%) reported here was more uniformly distributed down the profile than that reported for the global ecosystems (50%) [5] and under forests in China (54%) [26], possibly due to high OM inputs from both above-and belowground biomass, slow decomposition and stabilization of SOC by andic and spodic properties of non-volcanic andosolic soils [37].In the case of agricultural lands, the vertical distribution of SOC density (34%) is comparable to those reported for the cropland in China (37%) [26] but more homogeneously distributed compared to that under global ecosystems (41%) [5].As discussed above, the discrepancy may be due to relatively low SOC density in the upper layers, limited tillage under orchards and deep incorporation of farmyard manure and crop residues under dry land and paddy land.

Altitudinal Effects on SOC Density and Its Vertical Distribution
The clear association of LULC with both SOC and altitude means that there is a systematic variation in mean SOC density with altitude.The mean SOC density under different altitudinal zones was relatively large under 4000-5520 m (4.8-13.2kg•m −2 for first three depths) and 3500-4000 m zones (3.4-11.7 kg•m −2 ) compared to 3000-3500 m (2.2-11.6 kg•m −2 ), 2500-3000 m (1.7-7.0 kg•m −2 ) and 1769-2500 m zones (2.0-7.1 kg•m −2 ) (Table 6).The results from the mean separation tests showed significant differences (p < 0.05) in mean SOC density values between the top three and the last two altitudinal zones at 0-20 cm depth.Similarly, the values at 3500-4000 m and 3000-3500 m zones were significantly different (p < 0.05) from the last two zones at 20-40 cm depth.In the case of 40-60 cm depth, only 3500-4000 m zone was significantly different (p < 0.05) from the last two zones.The increase in SOC density with altitudinal zone was largely due to increase in SOC concentration with altitude resulting from higher OM inputs from above-and belowground biomass, slow decomposition due to low temperature [10] and more translocation of OC into deeper layers.The latter process was probably accentuated by increasing precipitation in the upper altitudinal zones (Table 1).The presence of non-volcanic andosolic soils in the upper altitudinal zones (3200-4000 m) might have stabilized SOC by forming complexes probably promoted by the soil andic and spodic properties [37].Previous studies have shown that interactions between SOC and soil minerals form organo-mineral complexes which in turn control the SOC dynamics and storage in the soil [40,41].
Conversely, the cumulative mean SOC density for the upper meter profile under different altitudinal zones decreased in the order of 3500-4000 m (33.31 kg•m −2 ) > 3000-3500 m (26.63 kg•m −2 ) > 2500-3000 m (17.17 kg•m −2 ) > 1769-2500 m zones (17.38 kg•m −2 ) with the 3500-4000 m zone significantly different (p < 0.05) from the last two zones.As discussed above, the increase in SOC density with increasing altitude of the altitudinal zone is largely due to increase in SOC concentration [13,40].In relative terms, the 3500-4000 m zone stores approximately 15.9 kg•m −2 (92%), 16.1 kg•m −2 (94%) and 6.7 kg•m −2 (25%) more SOC density in the upper meter profile than the 1769-2500 m, 2500-3000 m and 3000-3500 m zones, respectively.Our results were generally in line with previously reported mean SOC density values from the western Himalayas, India [42].A similar trend has also been reported by other studies [40,43].It was also found that the altitude factor affects the vertical distribution of SOC density in the profiles by influencing climatic factors like temperature and precipitation.A strong correlation of SOC with mean annual temperature and precipitation, and clay content has been reported [5].When altitude was plotted against SOC density, there was a positive correlation (r) of over 0.55 for the first three depths.The mean SOC density values under different altitudinal zones ranged from 7.1-13.2kg•m −2 (0-20 cm), 4.3-8.9kg•m −2 (20-40 cm), 2.1-5.0 kg•m −2 (40-60 cm), 1.9-3.0kg•m −2 (60-80 cm) and 1.7-3.4kg•m −2 (80-100 cm) (Table 6).Except for the 4000-5520 m zone, the mean SOC density values at 0-20 cm and 20-40 cm depths were significantly different (p < 0.05) from each other and also from the values in the last three depths.The other exception is that the mean SOC density values at 0-20 cm and 20-40 cm depths were also not significantly different under 3500-4000 m.The insignificant differences in SOC density values between 0-20 cm and 20-40 cm depths under 3500-4000 m and 4000-5520 m indicate that SOC density is more homogeneously distributed in the profiles compared to other altitudinal zones.As expected, the mean SOC density values in the last three depths were not significantly different and this suggests a decreasing effect of altitude on SOC density with depth.
The proportion of mean SOC density at 0-20 cm (relative to the first meter) under different altitudinal zones was highest under 3000-3500 m zone (43%) followed by 1769-2500 m and 2500-3000 m zones (41% each) and 3500-4000 m zone (35%) (Figure 3).The results indicate that the vertical distribution of SOC density under 3500-4000 m zone was more homogeneous than other zones possibly due to more dense forest cover in the upper altitudinal zones and presence of non-volcanic andosolic soils (3200-4000 m) with moderately high SOC content throughout the profile [30,37].In addition, slow decomposition and accumulation of SOC [10], more leaching and translocation of OC into deeper layers, and formation of organo-mineral complexes [40,44] such as with andic and spodic properties of non-volcanic andosolic soils [37] might have been the other contributing factors.

Impacts of Slope Aspect on SOC Density and Its Vertical Distribution
Slope aspect determines the micro-climate of hill/mountain slopes by influencing the solar radiation and evapotranspiration, and as such affects soil moisture, air temperature, vegetation cover and microbial activities.Aspect direction therefore impacts the SOC dynamics by its influence on the amount of OM inputs and subsequent decomposition in the soil.The mean SOC density values were variable in their ranges under different aspect directions, viz: 2.8-10.6 kg•m −2 (north-facing slopes), 2.6-9.0 kg•m −2 (east-facing slopes), 1.9-8.5 kg•m −2 (south-facing slopes) and 1.5-8.5 kg•m −2 (west-facing slopes) (Table 7).At all depths, the northern aspect had relatively large mean SOC density than other aspect directions.However, none of the aspect directions was significantly different (p < 0.05) at all depths.The cumulative mean SOC density in the top meter profile also depicted a similar trend with larger SOC density on the north-facing slopes (27.7 kg•m −2 ) than east (23.6 kg•m −2 ), south (21.4 kg•m −2 ) and west-facing slopes (17.1 kg•m −2 ).The larger SOC density on the north-facing slopes could be ascribed to rich biodiversity [19,45], high OM inputs [19,21], faunal abundance [19] and more soil moisture [19,20] on the northern aspect compared to other aspect directions.Our results are comparable with the reported mean SOC density of 21 kg•m −2 under north and 16 kg•m −2 under south-facing slopes (upper 1.1 m depth) in the subalpine range of the Italian Alps [20] and 4 kg•m −2 (south-western) and 18 kg•m −2 (north-western) in the first 30 cm depth in the temperate conditions of the Garhwal Himalayas, India [46].(1.5-2.8 kg•m −2 ) (Table 7).The mean SOC density values at 0-20 cm and 20-40 cm depths were significantly different (p < 0.05) from each other and from the values of the last three depths under all aspect directions.The exception is that the mean SOC density value at 20-40 cm depth was not significantly different from values at 0-20 cm and 40-60 cm depths under north-facing slopes.This suggests that SOC density on the northern aspect is more uniformly distributed down the profile compared to other aspect directions.On the other hand, the mean SOC density values in the last three depths were not significantly different (p < 0.05) except for values between 40-60 cm and 80-100 cm depths on the south-facing slopes.The relatively large SOC density in the top layer under all aspect directions is attributable to higher OM inputs from the aboveground biomass and more microbial abundance compared to the lower depths.Conversely, the proportion of mean SOC density in the upper 20 cm (relative to the cumulative mean SOC density in the upper meter) was highest under western aspect (49%) followed by southern (40%) and northern and eastern aspect directions (38% each) (Figure 4).The results indicate more homogeneous vertical distribution of SOC density (relative to the SOC density in the first 20 cm depth) on the north and east-facing slopes compared to the other two aspect directions.The larger and more uniform distribution of SOC density on the northern aspect could be largely attributed to rich biodiversity [19,21] and high OM inputs from above-and belowground biomass [19].Other factors include more soil moisture and leaching [20] and high clay content [47] for organo-mineral complexation [40,44] compared to other aspect directions.

Implication of Results
Change in LULC poses a threat to SOC storage in montane ecosystems.Soil erosion and land degradation would be inevitable as a result of LULC change particularly if natural vegetation is cleared and converted to agricultural lands.The results in this study show that if forest were to be converted to other LULC types, the first meter of the soil profile would release approximately 7.5 kg•m −2 , 7.0 kg•m −2 and 4.3 kg•m −2 SOC under agriculture, grasslands and shrublands, respectively.Since encroachment into forest land for agriculture is slowly becoming an issue in Bhutan and if conversion were to happen to a varying degree, this might trigger a rapid release of SOC from forest soils into the atmosphere.As such, appropriate land management practices and land use policies should be formulated to protect and conserve the current LULC types in order to ensure the natural forests remain a major C sink and provide continued ecosystem services.Moreover, the higher altitudinal zones with relatively large SOC storage are also susceptible to both anthropogenic influence and climate change.Fortuitously, climate change may lead to change in LULC type (e.g., shift in vegetation cover) along the altitudinal gradient; this might have a significant impact on biodiversity and SOC dynamics.In this regard, appropriate adaptation measures need to be considered to increase the resilience against climate change impacts.This study also shows marked variations in SOC density under different aspect directions with relatively large SOC storage on the north-facing slopes compared to other aspect directions.Given that LULC type, altitude and slope aspect have some significant impacts on the vertical distribution of SOC density in the profiles, incorporation of their effects in biogeochemical modelling would be useful for better understanding the C dynamics.

Conclusions
This study has shown some significant effects of LULC type, altitude zone and aspect direction on SOC density and its vertical distribution in the soil profiles.In the top meter, forests, 3500-4000 m zone and north-facing slopes stored the largest SOC density compared to other LULC types, altitudinal zones and aspect directions.The vertical distribution of SOC density under different LULC types was more homogeneous under agricultural lands than under forests, shrublands and grasslands.However, when forests, shrublands and grasslands were compared, forests had relatively large and more uniform distribution of SOC density down the profiles.Conversely, both 3500-4000 m zone and north-facing slopes had larger and more homogeneous vertical distribution of SOC density in the profiles compared to other altitudinal zones and aspect directions.Our results suggest that vegetation and soil moisture as influenced by altitudinal variation of precipitation are the main factors controlling the distribution of SOC density in the soil profiles.This has contributed to our understanding of the conditions favorable to most vigorous vegetative growth and litter production leading to deeper incorporation of OC in the soil.Results of this study have provided some insights of the potential losses of ecosystem services currently provided by the Himalayan forests, assuming its absence if they were to be cleared for other purposes.

Figure 1 .
Figure 1.Location of the study area in western Bhutan depicting sampling sites.

Figure 2 .
Figure 2. Proportion of SOC density (kg•m −2 ) stored at different depth intervals under different LULC types.Data are in percentages and are plotted at the midpoint of each depth interval.

Figure 3 .
Figure 3. Proportion of SOC density (kg•m −2 ) stored at different depths under different altitudinal zones.Data are in percentages and are plotted at the midpoint of each depth interval.

Figure 4 .
Figure 4. Proportion of SOC density (kg•m −2 ) stored at different depth intervals under different aspect directions.Data are in percentage and are plotted at the midpoint of each depth interval.

Table 1 .
Study site information under different altitudinal zones.

Table 2 .
Soil pH and bulk density (BD) (means ± standard error mean) at different depths under different LULC types.

Table 3 .
Effects of fixed terms (LULC, altitudinal zone and aspect direction) and their associated interactions on SOC density (kg•m −2 ) at each depth interval.

Table 4 .
Effects of fixed terms (LULC, altitudinal zone, aspect direction and depth) and their associated interactions on SOC density (kg•m −2 ) under each LULC type, altitudinal zone and aspect direction.

Table 5 .
may have been a contributing factor.Our results are comparable with the reported cumulative mean SOC density in the upper meter of 47 kg•m −2 under dense temperate forest and 22 kg•m −2 under open temperate forest in the Himalayas [38]; 14 kg•m −2 under evergreen needle forest, 12 kg•m −2 under deciduous broadleaf, 10 kg•m −2 under closed shrublands, 17 kg•m −2 under grasslands and 9 kg•m −2 under cropland [26]; and 17 kg•m −2 , 15 kg•m −2 and 12 kg•m −2 under temperate deciduous, evergreen and grasslands, respectively [5].Vertical distribution of mean SOC density values under different LULC types.
* Values in the columns suffixed with different superscript letters are significantly different (p < 0.05) from each other.

Table 6 .
Vertical distribution of mean SOC density values under different altitudinal zones.

Table 7 .
Vertical distribution of mean SOC density values under different aspect directions.
* Values in the columns suffixed with different superscript letters are significantly different (p < 0.05) from each other.