*Article* **Effects of Topography and Social Position on the Solar Radiation of Individual Trees on a Hillslope in Northwest China**

**Jiamei Li, Pengtao Yu \*, Yanfang Wan, Yanhui Wang, Bingbing Liu and Yipeng Yu**

Key Laboratory of Forest Ecology and Environment of National Forestry and Grassland Administration, Liupan Mountains Forest Ecosystems National Positioning Observation and Research Station, Ecology and Nature Conservation Institute, Chinese Academy of Forestry, Beijing 100091, China; lijiam@caf.ac.cn (J.L.); wanyf@caf.ac.cn (Y.W.); wangyh@caf.ac.cn (Y.W.); bingbing-l@caf.ac.cn (B.L.); yuyp@caf.ac.cn (Y.Y.) **\*** Correspondence: yupt@caf.ac.cn; Tel.: +86-10-62889562

**Abstract:** Solar radiation is a key factor influencing the photosynthesis and transpiration of trees. In mountainous regions, solar radiation income exhibits strong spatial heterogeneity due to topographical variations and the structural complexity of the forest. However, how the solar radiation income of individual trees in different social positions varies with slope position remains unclear. In this study, the daily solar radiation of the horizontal ground (*Rh*), different slope positions (i.e., at different locations on a hillslope, *Rs*) and individual trees with different social positions in the forest (*Ri*) were monitored from May to October in 2020 and 2021. The daily solar radiation income of a single hillslope (*Rf*) was applied to quantify the *Rs* response to the slope and aspect (i.e., slope effect) and the shade from the opposite mountain (i.e., shaded terrain effect). Our results showed that the *Rf* was 27.8% lower than *Rh* due to the slope effect of the sample slope. In the different slope positions, 2.7%–46.9% of solar radiation was lost due to the shaded terrain effect. A stronger limitation of *Rs* by the shaded terrain effect was detected on the bottom slope compared to that of the upper slope. The better the social position of an individual tree (i.e., tree dominance (*Dom*) and the distance between trees (*D*)), the more solar radiation it received, ranging from 22.4 to 95.3%. The dominant factor contributing to changes in *Ri* was slope position followed by *D* and *Dom* and, finally, *Rh*. These results provide an important basis for understanding the role of topography and tree social positions in solar radiation income in mountainous regions. Forest management measures should be varied with slope positions in mountainous regions, and forest density (i.e., distance between trees) should be considered as a key factor to optimize the forest functions.

**Keywords:** solar radiation; individual trees; topography; social position

#### Academic Editor: Francois Girard

**Citation:** Li, J.; Yu, P.; Wan, Y.; Wang, Y.; Liu, B.; Yu, Y. Effects of Topography and Social Position on the Solar Radiation of Individual Trees on a Hillslope in Northwest China. *Forests* **2023**, *14*, 561. https://doi.org/10.3390/

Received: 23 February 2023 Revised: 9 March 2023 Accepted: 9 March 2023 Published: 12 March 2023

f14030561

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **1. Introduction**

Solar radiation is one of the main environmental signals affecting plant biology, controlling multiple physiological responses [1,2], such as photosynthesis, transpiration and morphogenesis [3,4], with consequences for tree growth [5–8], which in turn can have an impact on forest productivity [9]. Solar radiation income often shows a high degree of spatial heterogeneity as a result of complex topographical variations and the spatial structure of trees [7,10–14], especially in mountain forests [15,16]. The accurate assessment of the effect of topography and the vertical and lateral position of trees in the forest (i.e., the niche of a tree in a forest is called social position in plant ecology) [17] on solar radiation income is essential for the measurement of solar radiation in mountain forests. Topography variations (such as slope gradient and aspect and the shadows cast by the surrounding topography) have a direct effect on solar radiation income [10,18–21]. At present, many studies of solar radiation in mountainous regions often focus on the difference in solar radiation income on various hillslopes [11,12,22]. For example, in the Northern Hemisphere, south-facing slopes

generally receive a greater amount of solar radiation [23]. Other studies have focused on topographic shading effects [24,25]. The shaded terrain effect by topography leads to a significant reduction in solar radiation income on hillslopes [16]. However, few studies have considered the differences in solar radiation income at different slope positions on hillslopes and whether the shaded terrain effect is consistent at different slope positions.

The solar radiation income of individual trees within forests on hillslopes is also affected by the tree's spatial structure (i.e., the spatial structure of the tree in both the vertical and horizontal orientation). The essence of the social position of individual trees can describe the spatial information of trees well, which is influenced by different environmental conditions for each tree (i.e., the relative size and distance between the reference tree and its neighbors) and, hence, individuals generate specific social position features [26,27]. The social position of trees determines the living space and available water and radiation resources. Trees with a dominant social position have a vast living space and plenty of water and radiation resources, causing less competition [17]. The degree of dominance and forest density can be used as a proxy indicator [9,28]. Some studies have suggested that the degree of tree dominance affects its ability to receive solar radiation [29], and dominant trees in a stand receive a greater amount of solar radiation [6,30]. Furthermore, decreasing the forest density increases the fraction of the absorbed solar radiation of trees [28,31,32]. However, in many studies, trees were simply divided into different dominance groups based on their relative height or diameter at breast height *(DBH)* in the forest [7,33–35], without considering the shading effects of the surrounding adjacent trees. Assessments of forest density also did not take into account the effect of the social position of individual trees on solar radiation income.

Forests are mainly located in mountainous regions [36]. The rugged mountainous terrain creates a major obstacle to forest management because of the difficult physical approach by men and machines [12,37]. Therefore, it is difficult to accurately quantify the effects of solar radiation income on mountain forests. In addition, forest management is currently focused on the stand scale, i.e., only the stand density is considered without considering the differences in the social location of trees, resulting in the forest management being invalid or only slightly valid [7,38].

To gain a deeper understanding of the effect of topography and social position on the solar radiation income of individual trees, we hypothesize that: (1) solar radiation income is similar at different slope positions; (2) solar radiation income is greater under a more dominant social position. Hence, we aimed to assess the effects of topography and social position on the solar radiation income in a mountainous region, which can help us understand the character of the solar radiation income of trees, so as to make targeted plans for forest management in mountainous regions.

#### **2. Materials and Methods**

#### *2.1. Study Site*

The study site was located in the small watershed of Xiangshuihe (*XSH*), Southern Liupan Mountains (106◦12 –106◦16 E, 35◦27 –35◦33 N) in Ningxia, Northwest China (Figure 1a). *XSH* has an elevation range of 2010 to 2942 m and a semi-humid climate, with a mean annual air temperature of 5.8 ◦C and an annual precipitation of 550.8 mm (1981 to 2010); 70% of precipitation occurs in the growing season from May to October. The soil is dominated by gray cinnamon with a sandy loam texture. *XSH* is a forested watershed with a forest cover of 82.4%, of which plantations account for 24.4%.

The plantations of larch *(Larix principis*-*rupprechtti*) are the main components of the plantations in *XSH*, accounting for 23.6% of the *XSH* area. *Larix principis*–*rupprechtii* is the main tree species used in afforestation in North and Northwest China. Such larch plantations play a very important role in producing timber and ecological services (e.g., soil protection, hydrological regulation and water supply) [39,40]. In larch plantations, the understory shrubs were scattered, covering only approximately 5% of the area. Herbs cover

approximately 40% of this area, with the dominant species being *Pteridium aquilinum* and *Carex hancokiana*.

**Figure 1.** Information about the locations of the study sample plots and sample trees. (**a**) Location of the Xiangshuihe (*XSH*) study watershed in NW China. (**b**) The location of plots (*US*, *MS*, *DS* and *BS* represent the upper, middle, down and bottom slope) along slope positions on the sample hillslope. (**c**) The detailed terrain of the study slope and opposite slope. (**d**) The valley cross-section.

#### *2.2. Study Slope and Plots*

#### 2.2.1. Study Slope

We selected a representative southeastern slope completely covered by an even-aged pure plantation of *Larix principis*-*rupprechtti*. The slope had a horizontal length of 425.1 m (corresponding to a slope length of 460.5 m), with a slope gradient of 22.6◦ and an elevation range of 2271 to 2480 m. The opposite mountain was the northwest slope with a horizontal length of 1052 m, a slope gradient of 23.5◦ and an elevation range of 2270 to 2731 m. The horizontal distance between the bottom slope of the opposite mountain and the bottom slope of the sample slope was 82 m (Figure 1b–d).

#### 2.2.2. Study Plots

In 2020, four plots were selected at different slope positions on the sample slope, vertically ranging from 2273 to 2471 m a.s.l. The upper slope (*US*), middle slope (*MS*), down slope (*DS*) and bottom slope (*BS*) were selected with elevations of 2471 m, 2397 m, 2293 m and 2273 m, respectively (Figure 1). These elevation points contain the whole range of the hillslope and the solar radiation received through these elevations can represent the difference amounts of solar radiation received on the whole slope. Also, since the shading of the opposite mountain was more obvious nearer to the bottom of the slope, measurement points were set up on both the down slope and bottom slope. The location information of the *Larix principis-ruprechtii* plantation plots on the study slope were shown in Table 1. The stand characters of each plot were shown in Table 2.


**Table 1.** The location information of *Larix principis*-*ruprechtii* plantation plots on the study slope.

The Δ*H* represents the elevation difference between the sample point and the opposite mountain top; the Δ*d* represents the horizontal distance between the sample point and the opposite mountain top.

**Table 2.** The stand characteristics of *Larix principis*-*ruprechtii* plantation plots on the study slope.


#### *2.3. Social Position of Trees in the Forest*

The social position of trees can represent the trees' vertical and lateral position in the forest. The trees in the forest were divided into four social position levels according to tree dominance (*Dom*) and distance between trees (*D*), as shown by Equations (1) and (2).

$$Dom = \frac{H\_i - \overline{H}}{\overline{H}} \tag{1}$$

where *Hi* (m) is the height of the sample tree *i* and *H* (m) is the average height of all the nearest neighbor trees around the sample tree [17]. The *Dom* represents the degree of dominance in the vertical direction of the sample tree from neighboring trees.

The distance among trees (*D*) is defined as the average horizontal distance between the sample tree and its neighboring trees, using Equation (2):

$$D\_i = \frac{1}{n} \sum\_{j=1}^{n} D\_{ij} \tag{2}$$

where *Dij* (m) represents the horizontal distance between the sample tree *i* and the neighbor tree *j*. The *Di* represents the degree of superiority or inferiority in the horizontal direction of the sample tree from neighboring trees.

In order to better describe the differences in the social position of trees within the forest, the sample plots at the down slope (*DS*) were used as an example to illustrate the characteristics of the social position of trees. There were 91 trees in the down slope plot, with a mean tree height (*H*) of 19.2 m, a mean diameter at breast height (*DBH*) of 20.63 cm and a mean crown diameter of 4.5 m (Table 2). All trees in the plot were divided into four social position levels according to tree dominance (*Dom*) and the distance between trees (*Di*). The vertical structure characteristics of the four social position levels are shown in Table 3.


**Table 3.** The vertical structure characteristics of *Larix principis*-*rupprechtii* plantation plot in the down slope (*DS*) (means ± standard deviations).

#### *2.4. Sample Tree Selection from the down Slope Plot (DS)*

In the down slope plot (*DS*), one to two representative trees from each social position level were selected as sample trees using a stratified random sampling method (i.e., sample trees that can represent the average of this level). Therefore, we selected a total of seven sample trees, and the *DBH* varied in the range of 15.50 to 25.92 cm, the tree height ranged from 17.2 to 21.8 m and the canopy diameter ranged from 3.5 to 5.2 m (shown in Table 4).

**Table 4.** The characteristics of seven sample trees in the down slope plot (*DS*).


#### *2.5. Solar Radiation Measurement*

2.5.1. Solar Radiation Data Conversion

The standard for measuring solar radiation utilizes the units of watts per meter squared (W m−2). However, pyranometers are both costly and limited in their ability to measure the solar radiation of individual trees in the slope field. With a lower cost, smaller size and more flexible installation, light meters measure luminous flux per unit area (illuminance), utilizing the units of lumens per meter squared or lux (lx). An effective conversion factor between watts per meter squared and lux would enable the use of light meters to evaluate the differences in the solar radiation in individual tree levels. Additionally, surveys of the literature have found no definitive and readily available "rule of thumb" conversion standard between solar radiation and illuminance data.

In this study, we converted the illuminance data (lx) into solar radiation data (W m<sup>−</sup>2) using Equation (3) (Figure 2). This equation was established based on data recorded from a pyranometer in an automatic weather station (Weatherhawk 232; WeatherHawk Inc., Logan, UT, USA) and a light meter (HOBO MX2202; Onset Computer Corp., Bourne, MA, USA). An automatic weather station and a light meter were installed at the same position (i.e., 1.3 m above the ground and the distance between them was 1.0 m) on open and flat horizontal ground approximately 100 m from the sample slope. The data concerning total solar radiation arriving on horizontal ground were collected every 5 min by an automatic weather station (Weatherhawk 232; WeatherHawk Inc., Logan, UT, USA) from May to October in 2020. Radiation data were collected every 5 min by a light meter (HOBO MX2202; Onset Computer Corp., Bourne, MA, USA) from May to October in 2020.

$$y = 0.009x - 1.360\tag{3}$$

where *y* (W m<sup>−</sup>2) represents the solar radiation and *x* represents the illuminance (lx).

**Figure 2.** Relationship between solar radiation and illuminance data.

#### 2.5.2. Solar Radiation of Individual Trees

The daily solar radiation income of individual trees (*Ri*, MJ m−<sup>2</sup> day<sup>−</sup>1) was measured with a light meter (HOBO MX2202; Onset Computer Corp., Bourne, MA, USA), which was installed at the top of the canopy of each of the seven sample trees. The sensor height of the sample trees is shown in the "Tree Height" column in Table 4. The data concerning the radiation arriving at the top of the canopy were collected every 5 min by the light meter (HOBO MX2202; Onset Computer Corp., Bourne, MA, USA) from July to October 2020 and May to October 2021.

#### 2.5.3. Solar Radiation over the Forest Canopy

The daily solar radiation income over the forest canopy (i.e., outside the forest) on the sample hillslope (*Rs*, MJ m−<sup>2</sup> day<sup>−</sup>1) was measured in each of the four slope plots (i.e., *US*, *MS*, *DS* and *BS*) with four light meter loggers (HOBO MX2202; Onset Computer Corp., Bourne, MA, USA) at a height of 1.3 m above the ground. The data concerning the radiation arriving at different slope positions were collected every 5 min by the light meter (HOBO MX2202; Onset Computer Corp., Bourne, MA, USA) from July to October 2020.

#### 2.5.4. Solar Radiation at the Horizontal Ground Measurement

The daily solar radiation income of the horizontal ground (*Rh*, MJ m−<sup>2</sup> day−1) was measured using the pyranometer of an automatic weather station (WeatherHawk 232, WeatherHawk Inc., USA), installed on open and flat horizontal ground approximately 100 m from the sample slope and at a height of 1.3 m above the ground [41]. The data concerning the total solar radiation arriving on horizontal ground were collected every 5 min by an automatic weather station (Weatherhawk 232; WeatherHawk Inc., Logan, UT, USA) from May to October in 2020 and 2021.

#### *2.6. The Calculation of Solar Radiation on a Single Hillslope*

The daily solar radiation income of the single hillslope (*Rf*) was introduced as a reference parameter (that is, the solar radiation of the slope is not affected by the surrounding mountains) and was calculated using the Fu (1958) model [10,42] as follows:

$$R\_f = (\mu \sin \delta + \theta \cos \delta \cos \omega - \sin \beta \sin \alpha \cos \delta \sin \omega) R\_h \tag{4}$$

$$\begin{aligned} \mu &= \sin\varrho \cos a + \cos\varrho \sin a \cos \beta \\ \theta &= \cos a \cos \varrho - \sin \varrho \sin a \cos \beta \end{aligned} \tag{5}$$

where *Rf* (MJ m−<sup>2</sup> day<sup>−</sup>1) represents the daily solar radiation income of the single hillslope, *β* ( ◦) represents the slope direction (expressed as the azimuth clockwise from north to east) and the slope gradient *α* ( ◦); *δ* ( ◦) is the declination; *ω* ( ◦) is the hour angle; and *ϕ* ( ◦) is

the latitude of the measurement point. The *u* and *v* are the composite parameters of the slope gradient *α* ( ◦), the slope direction *β* ( ◦) and the latitude of the measurement point *ϕ* ( ◦) using Equation (5).

#### *2.7. Relative Contribution of* Δ*H/*Δ*d, R, Dom and D to Solar Radiation Received by an Individual Tree*

The relative independent contributions of Δ*H*/Δ*d* (the ratio of the height difference to the horizontal distance between the sample point and the top of the opposite mountain, i.e., different slope positions, Δ represents the difference value), *Rh*, *Dom* and *D* to *Ri* were quantified using the coupled *Ri* model, and this model was constructed in the Results Section 3.4. The method for determining the relative contribution of each factor to *Ri* involved a comparison of the difference in the response variable (*Ri* in this study) between the simulation with only one varying factor (e.g., Δ*H*/Δ*d*, *Rh*, *Dom* and *D*) and the reference *Ri*.

The relative independent contributions of Δ*H*/Δ*d*, *Rh*, *Dom* and *D* to *Ri* were calculated using Equations (6) and (7) [39]:

$$\mathcal{C}\_{k} = \frac{\Delta R\_{i}(k)}{R\_{i}(reference)} \times 100\tag{6}$$

$$
\Delta \mathcal{R}\_i(k) = \mathcal{R}\_i(k) - \mathcal{R}\_i(reference) \tag{7}
$$

where *Ck* (%) is the independent contribution rate of each factor *k* (Δ*H*/Δ*d*, *Rh*, *Dom* and *D*) to the solar radiation of individual trees (*Ri*) compared to reference *Ri*. The reference *Ri* in this study was calculated using the developed *Ri* model and the value of Δ*H*/Δ*d* (0.0, i.e., the sample point has no shaded terrain effect), the long-term means of the *Rh* (9.32 MJ m−<sup>2</sup> day−1, 2013–2021) and the means of *Dom* (0.001) and *D* (2.01 m) from all the sample trees (835 trees) on the hillslope. This reference *Ri* value could represent the long-term mean environmental conditions and social position in this study hillslope. *Ri (k)* was calculated by inputting the measured value of factor *k* and the reference values of other factors into the *Ri* model.

#### *2.8. Statistical Analysis*

To clarify how *Ri* responds to each single factor (*Dom* and *D*) and the appropriate function type between *Ri* and *Dom* and *D*, the upper boundary line method [43], which can present the real relation between *Ri* and a single factor by eliminating the interferences of other factors, was used to determine the *Ri*–*Dom* and *Ri*–*D* relation. The upper boundary line was widely used in the establishment of compound models with independent factors as indicators [39].

Statistics were analyzed by the Statistical Product and Service Solutions (SPSS), version 19.0 (IBM Inc., Chicago, IL, USA) and presented as the mean ± standard deviation (SD). The differences in the solar radiation income at the varying slope positions and social positions was analyzed by one-way analysis of variance (ANOVA) with a significance level of *p* < 0.05, and the letters a and b indicate significant differences [41]. The performance of the model was assessed using the coefficient of determination (*R*2) [39].

#### **3. Results**

#### *3.1. Variation of Solar Radiation at Different Slope Positions*

The mean daily solar radiation income at different slope positions (*Rs*, slope gradient of 27.8◦, slope aspect of 135◦) was less than that of the horizontal ground (*Rh*) (Figure 3). The mean *Rs* at the upper slope (*US*) and the middle slope (*MS*) was slightly lower than the *Rh* (7.03 MJ m−<sup>2</sup> day−1) by 1.20 MJ m−<sup>2</sup> day−<sup>1</sup> and 1.60 MJ m−<sup>2</sup> day−1, respectively. However, *Rs* was even lower on the down slope (*DS*) and the bottom slope (*BS*), accounting for only 49.9% and 46.9% of *Rh*.

**Figure 3.** The differences in solar radiation income at various slope positions. (*Rh* represents the daily solar radiation income of horizontal ground; *US* represents the upper slope; *MS* represents the middle slope; *DS* represents the down slope; *BS* represents the bottom slope). 1.5 times IQR (interquartile range) was used to determine the anomalous values. Groups marked by different letters significantly differ from each other (i.e., a is significantly greater than b) (*p* ≤ 0.05).

Daily solar radiation income of the single hillslope (*Rf*) was significantly lower than the *Rh*, accounting for 72.2% of the *Rh* (Figure 4a). In addition, the *Rs* was always less than the *Rf* at the four slope positions by 2.7%–46.9%. The response trend of the *Rs* to the *Rf* could be expressed by a saturated logarithmic function, and the threshold of *Rs* decreased gradually with decreasing slope position (Figure 4b). When *Rf* was 8.46 MJ m−<sup>2</sup> day−1, the difference in solar radiation income between *US* and *MS* was not significant, and the corresponding thresholds were 8.23 MJ m−<sup>2</sup> day−<sup>1</sup> and 8.10 MJ m−<sup>2</sup> day−1, respectively. However, the difference in *Rs* was more pronounced at the *DS* and *BS* than at the *US*, and the corresponding threshold was 5.58 MJ m−<sup>2</sup> day−<sup>1</sup> and 4.49 MJ m−<sup>2</sup> day<sup>−</sup>1, respectively.

**Figure 4.** Effect of topographic change on solar radiation income on the hillslope. (**a**) The relationship between *Rf* and *Rh*. *Rh* represents the daily solar radiation income of the horizontal ground; *Rf* represents the daily solar radiation income of the single hillslope. (**b**) The relationship between *Rs* and *Rf* in different slope positions. *Rs* represents the daily solar radiation income outside the forest on the sample slope. The solid black line in (**b**) represents *Rf* = *Ri*. (**c**) The relationship between slope position and shade terrain effect. (*Rf*—*Rs*) represents the shaded terrain effect of the opposite mountain; Δ*H*/Δ*d* represents the ratio of the height difference to the horizontal distance between the sample point and the top of the opposite mountain (i.e., the variation of the slope position), the Δ*H* value represents the elevation difference between the sample point and the opposite mountain top; the Δ*d* value represents the horizontal distance between the sample point and the opposite mountain top; upper slope (*US*), middle slope (*MS*), down slope (*DS*) and bottom slope (*BS*). The (Δ*H*/Δ*d* < 0.26) was the threshold indicating that the sample point has no shaded terrain effect, that is, *Rf*—*Rs* < 0.5 MJ m−<sup>2</sup> day<sup>−</sup>1.

The response trend of *Rf*—*Rs* to Δ*H*/Δ*d* (i.e., the ratio of the height difference to the horizontal distance between the sample point and the top of the opposite mountain) could be expressed by a power function, and *Rf*—*Rs* first increased slowly with increasing Δ*H*/Δ*d* until 0.26 and then increased rapidly (Figure 4c). This indicated that the sample plot was not affected by the shaded terrain effect of the opposite mountain if the *Rf*—*Rs* was less than 0.5 MJ m−<sup>2</sup> day−1. In addition, the slope effect was always greater than the shaded terrain effect at *MS* and *US*, while the shaded terrain effect was more important than the slope effect at *DS* and *BS* (Figure 5).

**Figure 5.** Comparison of the slope effect and the shaded terrain effect at different slope positions. *Rslope* represents the slope effect (i.e., *Rh*—*Rf*).

#### *3.2. Effect of Tree Social Position on Its Solar Radiation*

The better the social position of an individual tree, the more solar radiation (*Ri*) it received (Figure 6). *Ri* at the dominant level was the largest, followed by *Ri* at the codominant level and intermediate level and, finally, the suppressed level. The mean *Ri* of the dominant level (4.82 MJ m−<sup>2</sup> day−1) was 1.1 times, 2.0 times and 2.4 times more than the codominant level, the intermediate level and the suppressed level, respectively.

**Figure 6.** The differences in the solar radiation income of individual trees among social position levels in the down slope plot. Groups marked by different letters significantly differ from each other (i.e., *a* is significantly greater than *b*) (*p* ≤ 0.05).

The *Ri* increased linearly with increasing daily solar radiation income from clearings outside of the forest on the sample slope (*Rs*) for the four social position levels, ranging from 22.4% to 95.3% (Figure 7a). The relationship between *Ri* at the dominant level and *Rs* was more sensitive than that of any other level. For example, the highest coefficient of determination (*R*2) of the dominant levels was 0.978, which was 1.0, 1.1 and 1.5 times higher than the codominant, intermediate and suppressed levels, respectively. The *UBLs* in Figure 7b show that *Ri* gradually increased with the increasing dominance of trees (*Dom*) following a linear function, namely, *Ri* = *e* × *Dom* + *f*. *Ri* first increased rapidly with increasing distance between the trees (*D*) until 4.5 m and then gradually stabilized around the maximum value as *D* continued to increase (Figure 7c). The relation between *Ri* and *D* could be expressed by a saturated logarithmic function, namely, *Ri* = *g* × ln(*D*) + *h*.

**Figure 7.** Effect of social positions on the solar radiation income of individual trees. (**a**) Response of *Ri* to *Rs*. (*Ri* represents the daily solar radiation income of individual trees; *Rs* represents the daily solar radiation income outside the forest of the sample slope). (**b**) Response of *Ri* to *Dom* (*Dom* represents the social position levels of the sample trees). (**c**) Response of *Ri* to *D*. (*D* represents the mean distance between the sample tree and neighboring trees). The upper boundary line (*UBL*) method [44] was used to determine the response function types of *Ri* to *Dom* and *D*.

#### *3.3. The Construction of the Tree Solar Radiation Model and Its Validation*

To understand the effect of topography and social position on the solar radiation income of individual trees (*Ri*), we built a solar radiation model of individual trees, which combined with the results of Sections 3.1 and 3.2. The *Ri* model was given as follows:

$$R\_i = R\_h \times f(topo\,\text{graph}) \times f(social\,\, position) \tag{8}$$

where *Ri* (MJ m−<sup>2</sup> day−1) represents the daily solar radiation income on the sample tree *i*, *Rh* (MJ m−<sup>2</sup> day−1) represents the daily solar radiation income of the horizontal ground, and *f* (*topography*) and *f* (*social position*) are the response functions of *Ri* to topography and social position factors, respectively. The response function of *Ri* to topography was obtained by the mathematical substitution method. The response functions of *Ri* to social position (i.e., dominance and distance) were obtained by concatenated multiplication since the solar radiation of individual trees showed a positive correlation with both dominance and the distance between trees (Figure 7b,c).

In addition, combining the methods in Section 2.6, the *Ri* model can be estimated by Equation (9):

$$R\_i = \left( a \times \left( (\mu \sin \delta + \theta \cos \delta \cos \omega - \sin \beta \sin \omega \cos \delta \sin \omega) R\_h - b \times \left( \frac{\Delta H}{\Delta d} \right)^{\varepsilon} \right) + d \right) \times \left( \varepsilon \times \text{Dom} + f \right) \times \left( g \times \ln(D) + h \right) \tag{9}$$

where *Ri* (MJ m−<sup>2</sup> day−1) represents the daily solar radiation income of individual trees, *β* ( ◦) represents the slope direction (expressed as the azimuth clockwise from north to east) and the slope gradient *α* ( ◦); *δ* ( ◦) is the declination; *w* ( ◦) is the hour angle; and *ψ* ( ◦) is the latitude of the measurement point. *u* and *v* are the composite parameters of the slope gradient *α* ( ◦), the slope direction *β* ( ◦) and the latitude of the measurement point *ψ* ( ◦) using Equation (5). *Rh* represents the daily solar radiation income of the horizontal ground, Δ*H*/Δ*d* represents the ratio of the height difference to the horizontal distance between the sample point and the top of the opposite mountain, *Dom* represents the social position levels of the sample trees, *D* represents the mean distance between the sample tree and neighboring trees and the *a*, *b*, *c*, *d*, *e*, *f*, *g* and *h* values represent the model parameters that need to be determined by the solar radiation data.

All parameters of the *Ri* model were newly fitted using observed data on the sample slope (slope 27.8◦, aspect 135◦) from July to October 2020, and the *Ri* model was further validated using the observed data from May to October 2021, Figure 8. The results indicated that the *Ri* model could still accurately estimate the varying *Ri* with a high *R*<sup>2</sup> value (0.83, n = 657). The simulated mean value (4.13 MJ m−<sup>2</sup> day<sup>−</sup>1) was 16.3% greater than the observed mean value (3.55 MJ m−<sup>2</sup> day<sup>−</sup>1). The *Ri* model of the sample slope was as follows:

$$R\_l = \left(0.296 \times \left(0.722 R\_{\rm h} - 107.617 \left(\frac{\Delta H}{\Delta d}\right)^{3.965}\right) + 0.950\right) \times \left(2.132 Dom + 2.032\right) \times \left(0.989 \ln(D) - 0.426\right) \text{ R}^2 = 0.83 \tag{10}$$

**Figure 8.** Comparison of the observed and simulated *Ri* values based on Equations (10)–(12) under the different weather conditions. (**a**) Comparison of the observed and simulated *Ri* values for all days. (**b**) Comparison of the observed and simulated *Ri* values for rain free days. (**b**) Comparison of the observed and simulated *Ri* values for all days. (**c**) Comparison of the observed and simulated *Ri* values for rainy days.

In addition, the *Ri* model accuracy of rain-free days was 0.01 greater than that of the weatherindependent model (Figure 8a), with *R*<sup>2</sup> = 0.84, and the accuracy of the *Ri* model fit was low for rainy days with *R*<sup>2</sup> = 0.61.

$$R\_i = \left(2.067 \times \left(0.722R\_h - 107.617\left(\frac{\Delta H}{\Delta d}\right)^{3.963}\right) - 1.909\right) \times \left(3.786Dom + 2.667\right) \times \left(0.789\ln(D) - 0.593\right) \text{ R}^2 = 0.84 \tag{11}$$

$$R\_i = (0.760 \times \left(0.722 R\_h - 107.617 \left(\frac{\Delta H}{\Delta d}\right)^{3.863}\right) + 0.314) \times (4.414 Dom + 0.825) \times (1.078 \ln(D) + 4.582) \text{ R}^2 = 0.61 \tag{12}$$

#### *3.4. Independent Contributions of* Δ*H/*Δ*d, Rh, Dom and D to Ri*

The relatively independent contribution of Δ*H*/Δ*d*, *Rh*, *Dom* and *D* to *Ri* compared to the long-term mean of reference *Ri* is shown in Figure 9. The Δ*H*/Δ*d*, *Rh* and *Dom* exerted slight negative effects on *Ri*, with contribution rates of −18.4%, −4.5% and −6.9%, respectively, while *D* had an obviously positive effect on *Ri*, with a contribution rate of 16.0%. This indicated that variations in slope position had a greater effect on *Ri* than social position and *Rh* (i.e., weather variations), and the dominant social position factor affecting the change in *Ri* was *D*.

**Figure 9.** Relative contribution rate of *Rh*, *Dom* and *D* to the daily solar radiation income of individual trees (*Ri*) compared to a reference *Ri*. The reference *Ri* was modeled by Equation (10) with the long-term means of *Rh*, the Δ*H*/Δ*d* value of 0 (Figure 5c) and the mean *Dom* and *D* of the 835 sample trees on the hillslope. (Δ*H*/Δ*d* represents the ratio of the height difference to the horizontal distance between the sample point and the top of the opposite mountain (i.e., the variation of slope position); *Rh* represents the daily solar radiation income of the horizontal ground; *Dom* represents the social position levels of sample trees; *D* represents the mean distance between trees and neighboring trees).

The response of *Ri* to *D* exhibited the same trend at different slope positions and different *Dom* levels, and an obviously higher rate of increase occurred under higher *Dom* levels (Figure 10). For example, given the Δ*H*/Δ*d* value of 0.1 and *Dom* value of 0.3, the effect of *D* on *Ri* only gradually decreased when *D* was greater than 6.0 m, while given the Δ*H*/Δ*d* value of 0.1 and *Dom* value of −0.3, *Ri* stabilized when *D* was greater than 3.0 m.

**Figure 10.** Variation in the simulated *Ri* with *Rh* at different *Dom* and slope position levels by the *Ri* model (Equation (10)). The value of *Rh* of 9.32 MJ m−<sup>2</sup> day−<sup>1</sup> used the long-term mean values from 2013 to 2021. Δ*H*/Δ*d* represents the ratio of the height difference to the horizontal distance between the sample point and the top of the opposite mountain. Δ*H*/Δ*d* = 0.26 was the threshold indicating that the sample point had no shaded terrain effect (Figure 4c).

There was a significant difference in the effect of *D* on *Ri* under different slope positions (i.e., different Δ*H*/Δ*d* values) (Figure 10). The difference in *Ri* at different *Dom* levels gradually decreased as Δ*H*/Δ*d* values increased, especially after Δ*H*/Δ*d* > 0.26. For example, given the Δ*H*/Δ*d* value of 0.1 and *Dom* value of 0.3, the *Ri* was above the multi-year means (9.32 MJ m−<sup>2</sup> day−1) if the *D* values were higher than 4.0 m. However, given the Δ*H*/Δ*d* value of 0.45, the *Ri* was always less than 9.32 MJ m−<sup>2</sup> day−<sup>1</sup> regardless of the value of *Dom*. This indicated that the higher the Δ*H*/Δ*d* value (i.e., the lower the slope position), the weaker the effect of the social position of trees (both *Dom* and *D*).

#### **4. Discussion**

#### *4.1. The Role of Topography on Solar Radiation: Slope Effect and Shaded Terrain Effect*

In this study, the solar radiation income on the sample slope (slope 27.8◦, aspect 135◦) was lower than that of the horizontal ground (Figure 3). Similar results were obtained in earlier studies [45,46]. This was mainly because of the slope effect (i.e., slope and aspect) and the shaded terrain effect caused by the opposite mountain [25,47]. Since the slope effect of solar radiation income was the same on a hillslope [42,48], the difference in solar radiation income at different slope positions was mainly caused by different shaded terrain effects. Therefore, we rejected the first hypothesis that the solar radiation income was similar at different slope positions.

The difference in the shaded terrain effect was due to the different blocking effects of topographic features (that is, the distance and elevation difference from the sample slope) [21,49,50]. Our research showed that the shaded terrain effect did not limit solar radiation income on the slope when above the middle slope (i.e., the ratio of the elevation difference to the distance between the sample point and the top of the opposite mountain was satisfied (Δ*H*/Δ*d*) < 0.26) (Figure 4). This was mainly because the solar radiation income above the mid-slope was not affected by mountain shading due to the difference in height at *US* and *MS*, with the opposite mountain being reduced and the distance between the opposite mountain being increased. This was consistent with the result indicating that the forest above the middle slope was not affected by shaded terrain [37].

#### *4.2. The Distance among Trees Was the Dominant Factor for Tree Solar Radiation*

The social position of trees (including dominance of trees and distance between trees) was unique for each tree in the forest [14]. Our results suggested that the higher the social position of the individual tree, the more solar radiation it received (Figure 6). This result was completely consistent with our second hypothesis. This was primarily because the trees located in the dominant canopy, which were less shaded by neighboring trees because of a relatively higher tree height [51], could receive more solar radiation [6,7]. Solar radiation income increased as the distance between trees increased (Figure 7), which was consistent with the result indicating that an increase in the distance between trees resulted in greater intakes of solar radiation [52].

However, the contribution of the distance between trees to *Ri* was more important than that of dominance (Figure 9) because the vertical shading effect among trees weakened as the distance between trees increased [53]. Moreover, the solar radiation income of individual trees increased slowly when the distance between trees was above 4.5 m (Figure 7c), indicating that these trees obtained sufficient solar radiation and the distance between trees was no longer a key factor limiting solar radiation income. This was consistent with the finding that competition among neighbors within5mreduced the growth rate due to limited light resources [7,54].

#### *4.3. Applications of the Newly Developed Ri Model*

The *Ri* model developed in this study exhibited a satisfactory performance in estimating the daily solar radiation income of individual trees. The simulated value was 16.3% higher than the observed value since the impacts of leaf variation on solar radiation income were not considered [44]. When simulating only rainy days, the model accuracy was lower at *R*<sup>2</sup> = 0.60 (Figure 8c). This was mostly because on rainy days, cloud cover and ground humidity (changing albedo) in turn affected solar radiation income [55,56]. The developed *Ri* model couples the effects of topography (including latitude, longitude, elevation, slope and aspect) and the social position of trees (both horizontal distance and vertical dominance), which can easily be obtained in different mountain regions. Foresters could choose suitable measurements according to their actual conditions. For example, horizontal ground solar radiation data can be obtained from ground weather stations, topographic data can be obtained from remote sensing monitoring and the social position of trees can be obtained through field surveys [57].

The *Ri* model can be used to explain tree growth and thereby regulate forest productivity. For example, a study on this sample slope showed that the cumulative seasonal growth of the down slope plot (1.08 mm) was significantly less than that of the upper slope plot (1.54 mm) because high solar radiation facilitates photosynthesis of organic matter for tree growth and increases air temperature [58]. In addition, on this sample slope, previous transpiration studies focused on the stand scale because of the lack of environmental characteristics of individual trees [37,59,60]. The *Ri* model in this study accurately predicted the solar radiation received by individual trees, which can be beneficial in explaining the differences in the transpiration of trees. In summary, a series of mountain forest management plans (e.g., thinning) can be made based on the topographic conditions and social position of trees using the developed *Ri* model to improve the forest structure and alleviate the adverse effects of shade on the solar radiation income of trees.

#### *4.4. Implications for Mountain Forest Management*

The solar radiation income of individual trees in forests depended on topographic and social positions. Therefore, different management measures should be adopted for trees at different slope positions and different social positions to promote the receipt of more solar radiation. Our research suggests that the social position of trees (especially the distance between trees) has relevance for mountain forest management.

Different forest management measures should be implemented at different slope positions. When Δ*H*/Δ*d* > 0.26, the solar radiation income of these slope positions was significantly affected by the shaded terrain effect of the opposite mountain; therefore, the forest density should be reduced by increasing the distance between trees so that residual dominant trees can obtain more solar radiation during the hours of sunshine. For example, given a Δ*H*/Δ*d* of 0.40 and a *Dom* of 0.3, the trees can receive a maximum *Ri* (8.79 MJ m−<sup>2</sup> day<sup>−</sup>1) at a distance between trees of 10.0 m (Figure 10).

However, when Δ*H*/Δ*d* < 0.26, the shaded terrain effect was almost negligible for the forest, and the shading of adjacent trees was the main factor affecting the solar radiation income of individual trees. The proportion of trees with different social positions can be adjusted by implementing thinning to reduce the light competition among trees. For example, given a Δ*H*/Δ*d* of 0.10, the daily *Ri* will be above 9.32 MJ m−<sup>2</sup> day−<sup>1</sup> when the *Dom* values are 0.2 and 0.3 with *D* values of 6.5 m and 4.0 m, respectively (Figure 10).

#### *4.5. Limitations of This Study*

The *Ri* model established in this study can calculate the solar radiation income of trees at any location on the slope. However, due to the experimental conditions, this study only monitored the actual solar radiation income of trees on a sample down slope, so the differences in the solar radiation income of trees in different slope positions could not be verified. Therefore, the reliability of the model estimation will be improved when the model parameters are fitted using the data from more mountainous regions [61].

#### **5. Conclusions**

This study focused on the effects of topography and social position on the solar radiation income of individual trees. The slope effect resulted in lower solar radiation income on the hillslope, which was 27.8% lower than that of the horizontal ground. In the different slope positions, 2.7% to 46.9% of solar radiation was lost due to the shaded terrain effect. Solar radiation income above the middle slope was basically unaffected by the shading of the opposite mountain. The more dominant the tree's social position, the more solar radiation it received, ranging from 22.4% to 95.3%. The dominant factor contributing to the changes in *Ri* was the distance between trees, and *Ri* was basically stable when the distance between trees was greater than 4.5 m. The *Ri* model covered the effect of topography and the social position of individual trees, which explained 83.0% of the observed *Ri* variation.

These results advanced our understanding of the role of topography and tree social positions of solar radiation income in mountainous regions. Therefore, different forest management strategies should be developed at different slope positions, and the distance between trees should be selected as the main management measure to facilitate the precise adjustment of stand structure design, ensuring that residual trees receive more solar radiation.

**Author Contributions:** Formal analysis, J.L.; investigation, J.L. and B.L.; writing—original draft, J.L.; writing—review and editing, J.L., P.Y., Y.W. (Yanfang Wan), Y.W. (Yanhui Wang) and Y.Y. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was financially supported by the National Natural Science Foundation of China (U21A2005, 42161144008, U20A2085) and the National Key R & D Program of China (2022YFF0801803, 2022YFF1300404).

**Data Availability Statement:** Not applicable.

**Acknowledgments:** We thank Liu Hong, Songping Yu, Zhonghui Zhang and Xiaonan Huang for their assistance in the field work.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

**Dmitry Pershin 1,2,\*, Natalia Malygina 2, Dmitry Chernykh 2, Roman Biryukov 2, Dmitry Zolotov <sup>2</sup> and Lilia Lubenets <sup>2</sup>**


**\*** Correspondence: pershindk@my.msu.ru

**Abstract:** The stable water isotopes in snow (primarily 18O and 2H) are widely used for tracing hydrological and ecological processes. However, isotopic signatures of snow can be significantly modified by topography and land cover. This study assesses spatial and temporal variability of the bulk snowpack isotopic composition (δ18O, δ2H, d-excess) between forested (pine and birch) and open areas in the West Siberian forest steppes. Isotopic samples were collected over the peak snow accumulation in 2017–2019. The snow isotopic composition within forested areas differed from open steppes, mainly in reducing d-excess (1.6‰ on average). We did not find a significant effect of canopy interception on snow enrichment in heavier isotopes. Snowpack in the pine forests was even lighter by 3.6‰ for δ2H compared to open areas, probably, due to low energy inputs and interception capacity. Additionally, snow depth significantly influenced the isotopic composition spatial variability. As snow depth increased, δ18O and δ2H values decreased due to conservation within the snowpack and less influence of sublimation and moisture exchange with the soil. However, this pattern was only evident in winters with below-average snow depth. Therefore, taking into account snow depth spatial and seasonal variability is advisable when applying the isotopic methods.

**Keywords:** river basin; forest; grassland; interception; wind redistribution; stable water isotopes

#### **1. Introduction**

Snow accumulation, storage, and melting dynamics affect multiple hydrological, ecological, and social processes in mountainous and high-latitude environments [1,2]. Tracking changes in snowpack accumulation and melt rates is challenging because the driving factors operate at multiple spatial and temporal scales [3–5]. Stable water isotopic composition of snow (primarily δ2H and δ18O) has become a valuable tool for investigating various snow hydrological and ecohydrological processes [6]. Implementations include snow contribution to groundwater recharge [7,8], streamflow generation during rainon-snow events [9,10], exploring vegetation water sources [11–14], etc. However, the application of isotopic methods is complicated by snow evolution processes that alter the isotopic signal [15–17].

Snow mass and energy balance is altered by complex processes such as sublimation, wind redistribution, forest canopy interception, melting, and metamorphism [18]. Since snow contains liquid, solid, and vapor phases, most of these processes are accompanied by phase transitions, changing the stable water isotopic composition [6,19,20]. The intensity of snow hydrological processes varies in space. Spatial factors affecting snow isotopic composition have included altitude [21–23], aspect [15,24], snow depth [21], and canopy interception [25–27]. Most of these factors affect snow sublimation fluxes, which leads to enrichment in heavier isotopes of the remaining snow cover [16,22,28,29].

Forest canopy interception affects the snow isotopic composition by increasing the sublimation of intercepted snow and making throughfall isotopically heavier [6]. The undercanopy snowpack in the north-western US was up to one-fourth smaller and isotopically

**Citation:** Pershin, D.; Malygina, N.; Chernykh, D.; Biryukov, R.; Zolotov, D.; Lubenets, L. Variability in Snowpack Isotopic Composition between Open and Forested Areas in the West Siberian Forest Steppe. *Forests* **2023**, *14*, 160. https:// doi.org/10.3390/f14010160

Academic Editors: Yanhui Wang, Karl-Heinz Feger and Lulu Zhang

Received: 23 November 2022 Revised: 12 January 2023 Accepted: 13 January 2023 Published: 16 January 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

heavier by roughly 2‰ in δ18O compared with the snowpack in the clear-cut area [26]. Studies in Switzerland showed that the isotope ratios were higher in the snowpack under forest canopy than in open grasslands (by 13.4 ‰ in δ2H and 2.3 ‰ in δ18O) [27]. However, a five-year study in the southwestern US has shown a much more significant influence of snowfall isotopic input and aspect on δ18O than canopy density [28]. Additionally, several mechanisms of changes in the snow isotopic composition in complex landscapes remain poorly understood, such as the effect of wind redistribution [6].

Most of the works cited above were performed in mountain forests in relatively humid regions of Europe and North America. Over Siberia, studies have shown that the variability of the contribution of many precipitation sources to the snow results in large isotopic variability [30,31]. At the same time, more detailed catchment-scale studies investigating changes in snow isotopic composition have not been conducted either in the boreal forest area or in the forest steppe.

Snow in continental semi-arid regions is the major water source for ground and soil water recharge and streamflow generation [32]. Considering the significant differences in the snow isotopic signal compared with rainfall, implementation of the isotopic methods for tracing streamflow formation, groundwater recharge, and plant water use seems promising in these regions. However, the mechanisms of the snowpack isotopic composition spatial variability and post-depositional fractionation remain poorly understood. In this work, we focused on changes in the stable water isotopic composition of snow (δ18O, δ2H, d-excess) between open and forested areas during peak snow accumulation over three years (March 2017–2019). The studies were conducted in the Kasmala River basin located in the forest steppe ecoregion in the south of Western Siberia. The basin structure consists of extensive arable lands, Scots pine forests, and small patches of deciduous forests. Basin landscape composition allowed us to study changes in the snow isotopic composition in open and forested areas and to evaluate factors influencing isotopic ratios considering the differences in topography and seasonal climate.

#### **2. Materials and Methods**

#### *2.1. Study Area*

We carried out the study in the 1768.7 km2 Kasmala river basin (53◦4' N, 82◦20' E) in the south of the West Siberian Plain (Figure 1). The Kasmala basin is a snow-dominated watershed in the headwaters of the Ob River.

The study area belongs to the West Siberian forest steppe ecoregion [33]. This region is also part of the West-Siberian grain belt, an important agricultural region in southern Siberia. The dominant land cover type is arable land (59.7%). A unique characteristic of this part of the West Siberian Plain is the long strips of Scotch pine (*Pinus sylvestris* L.) forests [34–36]. This forest type (covers about 12%) is characterized by dense pine stands (canopy density 60–80 %) with a small proportion of mixed deciduous vegetation. The average diameter at breast height is about 25–30 cm, and tree height is between 20–25 m. The forests occupy the sandy massifs within the extended ancient flow depressions. Upland slopes oriented towards the depressions have small slope angles (1–3◦) covered by arable land and steppe patches. Another type of forest presented here is the small patches of birch (*Betula pendula* Roth) and aspen (*Populus tremula* L.) stands (5.9%). The basin is naturally divided into three main parts: the northern (NP) and southern (SP) open steppe areas with deciduous forest patches and the central part, occupied mainly by pine forest (CP). Additionally, in the study, we considered arable land/open steppes, deciduous forests, and pine forests as three major land cover types.

**Figure 1.** The Kasmala basin in the south of Western Siberia with sampling locations and land cover composition. NP—northern part, CP—central part, SP—southern part.

Continental climate with strong seasonality in the study region varies from sub-humid in the northeast to semi-arid (steppe) in the southwest. The area has cold winters and hot summers with a mean annual temperature (1966–2020) of 2.4 ◦C at the Barnaul weather station (53 km to the nearest transect) [37]. The average annual precipitation is 427 mm. The mean winter (November-March) precipitation is about 125 mm. Wind drift is the most crucial factor in snow redistribution. The average daily wind speed in November-March is 2.7 m/s [37].

#### *2.2. Sampling and Laboratory Analysis*

We conducted three field campaigns around 10–15 March, 2017–2019 during the peak snow accumulation. Sampling was performed along eight transects that varied in length from 500 m to 2 km. The distance between measurements was 200 m in open areas and 100 m in forested areas. NP and SP included three transects each, and CP included two transects (Figure 1). Each transect differed in terrain characteristics (ruggedness, slopes, and aspects) and land cover composition (open areas, deciduous, and coniferous forests). GPS location of each sampling point was recorded. During the field campaign, the bulk snowpack samples were collected at each location using a 60 cm snow coring sampler VS-43. The samples were weighed to obtain bulk density and calculate snow water equivalent (SWE). Then, the entire samples were immediately placed into sealable high-density polyethylene bags to prevent further evaporation and fractionation effects. The snow depth was measured at the same points using a special snow ruler. If the snow depth was higher than 60 cm, the snow was measured successively from deeper layers. A total of 192 samples were taken over three sampling years.

Snow samples were brought to the laboratory in sealed bags on the day of sampling and melted at room temperature. Once the samples completely melted, they were filtered through 0.45-μm filters (Minisart NML Plus) into 2-mL glass vials. The isotopic composition (δ18O, δ2H) of all samples was analyzed through laser spectroscopy (PICARRO L2130-i (WS-CRDS). The measurement uncertainties were ±0.4‰ for <sup>δ</sup>18O and ±0.1‰ for <sup>δ</sup>2H. Water samples were calibrated against the international standards (V-SMOW, GRESP).

#### *2.3. Data Processing and Analysis*

We analyzed spatial and temporal variations in the bulk snowpack isotopic composition. In exploring interannual differences, we considered variations in temperature, precipitation, and snow accumulation relative to the interannual means from the Barnaul weather station [37]. The station is the closest to the study area (53 km to the nearest transect) and has the most consistent series of observations.

The δ18O and δ2H mean values among basin parts, land cover types, and sampling years were tested using non-parametric Kruskal–Wallis H test. Additionally, Wilcoxon test with the Bonferroni correction was applied to evaluate the differences between the group levels.

For assessing the site-specific covariation of hydrogen and oxygen stable isotope ratios, we created δ2H vs. δ18O plots (similar to local meteoric water lines) for the basin parts (NP, CP, and SP). Additionally, the second-order isotopic parameter deuterium excess was computed (d-excess = <sup>δ</sup>2H–8 ×δ18O) [38]. D-excess helps distinguish equilibrium and nonequilibrium processes through the differences from the global meteoric water line (GMWL; δ2H = δ18O + 10) [39]. Values less than 10‰ often indicate kinetic fractionation due to evaporation or sublimation [6].

In order to assess the interactions between isotopic composition, topography, and land cover, we created a set of multiple regression (MLR) models. Stepwise regression analysis was performed separately for δ18O and δ2H as dependent variables. However, the relationships behaved similarly, and we used δ2H for further analysis due to a smaller measurement error. All variables were first tested for normality of distribution, multicollinearity, and the presence of outliers. We also calculated Z-scores for the predictors and response variables of the final models. The selection of variables for each model was performed automatically by stepwise backward elimination according to the Akaike information criterion (AIC). Finally, each final model was checked for residual distribution and homoscedasticity.

Input variables were chosen based on their potential effect on the isotopic composition of a bulk snowpack. Snow depth and SWE were used as indicators of wind redistribution since wind is the main factor of depth and SWE spatial heterogeneity on the plains (no elevation gradient). Given that depth and SWE are correlated, they were included in the regression separately. The snow energy balance and, therefore, the processes of sublimation and melt depend on the local topography. The slope, aspect, general curvature, and Topographic Ruggedness Index (TRI) were calculated using the SRTM digital elevation model (DEM) [40]. The aspect was transformed with the cosine function, resulting in a value of −1 for north- and 1 for south-facing points. We used the land cover map based on the Landsat satellite images to calculate the forest ratio in the 200 m surrounding the sample points. This metric was chosen because forests in the study area affect snowpack accumulation not only via canopy interception but also through the obstruction of wind redistribution. The obstruction factor appears mainly in the deciduous forest patches within the steppes.

Statistical analysis was performed in R (http://www.r-project.org, accessed on 27 March 2021). Topographic variables were calculated using QGIS (www.qgis.org, accessed on 2 June 2021).

#### **3. Results**

#### *3.1. Interannual Differences*

Most seasonal climate parameters were close to average during winter seasons 2016/ 2017–2018/2019, excluding snow conditions. Mean winter temperatures did not differ much from the long-term mean (Table 1). Daily temperatures above 0 ◦C) were uncommon (for example, in December 2018). However, very low temperatures occurred more

frequently, such as in November 2016, January 2018, and February 2019 (Figure 2). Wind speeds in 2016/2017 were close to the mean (2.3 m/s), 2017/2018 and 2018/2019 had lower average wind speeds (1.7–1.8 m/s).

**Table 1.** Winter summary statistics with standard deviations at Barnaul weather station calculated from daily data [37].


**Figure 2.** Barnaul weather station air temperature, wind speed, snow depth for the studied winter seasons as well as the station period of record (1966–2020).

The strongest differences were observed in snow cover dynamics. The maximum snow depth values were above average in 2016/2017 (73 cm) and extremely below average in 2017/2018 (25 cm). March 2018 had the lowest snow depth during the station period of record. Depth and SWE values also varied on our transects during the study period. The most considerable difference occurred between 2016/2017 and 2017/2018, when depth and SWE were 57 % and 70 % lower, respectively (Table 2). Therefore, the studied winter seasons had significantly different snow mass conditions but were relatively similar in other seasonal climate parameters.

**Table 2.** Means and standard deviations of the bulk snowpacks SWE, depth, and isotopic composition, averaged over 2017–2019 and northern (NP), central (CP), and southern (SP) parts.


According to Kruskal–Wallis test results (Table 3), differences in mean isotopic composition between sampling years were significant. The median levels varied only around 1‰ in δ18O and 4‰ in δ2H during sampling years.



Variability within each sampling year was substantially higher (Figure 3). The isotopic variation had the highest rates in 2018 when the snowpack was extremely shallow. The range reached 6.3‰ in δ18O and 42.6‰ in δ2H. In 2017, the range was approximately half and amounted to 3.2‰ and 26.0‰ in δ18O and δ2H, respectively. The standard deviations (Table 2) of each sampling year also demonstrate differences in variability.

Differences in mean isotopic composition among the northern (NP), central (CP), and southern (SP) parts of the basin were also significant considering the whole study period (Table 3). Firstly, significant differences may indicate expected differences in the drivers of snow isotopic composition: primarily wind redistribution in open areas and canopy interception in forests. Looking at each year separately, the forested and open parts of the basin differed significantly only in δ2H (Figure 3). Additionally, Wilcoxon test showed (*p*-value < 0.05) that δ18O and δ2H varied significantly between NP and CP-SP.

**Figure 3.** Temporal variability of the isotopic composition (δ2H and δ18O) of bulk snowpack over the basin parts (NP, CP, SP) with results of Kruskal–Wallis tests within sampling years. Violin and Tukey outlier box plots show the median value (horizontal line within the box), the 1st and 3rd quartile (ends of the box), minimum/maximum values (whiskers) and distribution.

Differences among land cover types were significant for δ2H and d-excess mean values (Table 3). As in the case of the basin parts, δ2H and d-excess varied significantly between open steppes, deciduous forests, and pine forests (Wilcoxon test *p*-value < 0.05). Pine forests and deciduous forests were not statistically different.

During the entire study period, the snowpack isotopic composition was slightly heavier within the open parts compared with the forested. The δ18O and δ2H values in CP were on average lower than in open NP and SP, but only by 0.2 and 3.6‰, respectively. We expected higher δ18O and δ2H values in the forested part due to canopy interception. However, no direct effect was observed.

#### *3.2. Oxygen and Hydrogen Ratios*

Figure 4 gives an overview on δ2H vs. δ18O relationships among the sampling years and basin parts. All samples lay lower than the GMWL, indicating the influence of nonequilibrium processes (Figure 4). In all equations, slopes were far below 8, and intercept values did not exceed 10, corresponding to the similar parameters of GMWL.

Slope and intercept values of δ2H vs. δ18O regressions showed clear spatial patterns. In the open SP and NP slopes ranged from 7.90 to 6.12. The slope values were significantly lower in the forested part (5.58 to 6.38). Furthermore, the CP intercepts also reached very low values (up to –43.3). This aspect indicated significantly higher sublimation rates and other nonequilibrium processes in the forested part.

Slopes changed little over time, especially depending on the amount of snow each winter season. Only the intercept values changed significantly (by more than 10–20) within the open and forest parts. On the one hand, this may indicate the constancy of fractionation processes on-ground (stable slope), and on the other hand, the influence of snow formation conditions during each winter season (unstable intercept).

**Figure 4.** δ2H vs. δ18O plot over the basin parts (NP, CP, SP) and sampling years, including their regression lines (blue), global meteoric water line (gray), and snow depth gradient.

D-excess values also showed spatial patterns, remaining relatively stable interannually (Figure 5). CP's mean d-excess values (Table 2) were 1.6 ‰ lower than open SP and NP. Moreover, the d-excess fell below 10 (d-excess of GMWL) in all parts of the Kasmala basin. According to the Kruskal–Wallis test (Table 3), mean differences were significant across the basin parts during the study period. The major differences occurred between CP and NP/SP (Wilcoxon test *p*-value < 0.001). The contrast again referred to the high intensity of nonequilibrium processes in the forested part of the basin, which, however, did not lead to direct enrichment in heavier isotopes of the forest snowpack.

**Figure 5.** D-excess values over the basin parts (NP, CP, SP) and sampling years. Individual samples are shown using dots. Tukey outlier box plots show the median value (horizontal line within the box), the 1st and 3rd quartile (ends of the box), and the minimum/maximum values (whiskers).

#### *3.3. Influence of Topography and Land Cover Factors*

Only snow depth and TRI had a relatively stable effect on δ2H. The influence was evident in two of the three years. The other predictors had no significant influence, except for a feeble impact of Aspect in 2017 (Table 4).

**Table 4.** Adjusted R2, slopes, and intercepts for computed regression equations describing H2 vs. land cover and topography predictors relationships (ns = predictor variable not selected).


\* *p* < 0.1; \*\* *p* < 0.05; \*\*\* *p* < 0.01.

The effect of snow depth was closely related to the snow amount over the sampling years. The influence of depth was not significant in 2017 (high snow). However, in 2018 (shallow snow), the influence of depth became considerable and described almost 32% of the variability. The slope was negative, which means that δ2H values decrease with increasing snow depth. In moderate snow 2019, the influence of depth became weaker. The linear influence of snow depth is evident in the plot of partial residuals (Figure 6).

**Figure 6.** Partial residuals of each predictor for the sampling years selected by the stepwise multiple regression. The dashed lines represent linear relationships.

The TRI effect was most likely related to the wind redistribution of snow from uplands to small valleys. Higher TRI values mean higher ruggedness of the terrain. In the study area, higher TRI values typically belong to river valleys. The coefficients were also negative, indicating a decrease in δ2H values as we moved toward the valleys. TRI was a significant predictor in 2017 and 2019. These winters had above-average snow depths and quite a high wind intensity, especially in 2017.

SWE, as well as forest ratio, did not show any significant effect on the isotopic composition. SWE was included in the models instead of snow depth, but it was significant only in 2018, and the influence was weaker (adjusted R2 was about 0.15). A similar SWE value can be formed by increasing both depth and density, which are controlled by distinct processes from the isotopic composition perspective. The forest ratio also did not play a significant role. We attribute this effect to fundamentally different processes occurring in coniferous (canopy interception) and deciduous forests (obstacle of wind transport). However, the sampling points in these areas may have similar forest ratio values.

The d-excess also showed strong correlations with snow depth and SWE, confirming the isotopic composition-depth patterns (Figure 7). Snow depth and d-excess were positively correlated (Table 5) in 2019 (moderate snow) and 2018 (shallow snow). This relationship weakened in the high snow winter season (2017). The relationship was most stable in moderate 2019, while in 2018 several points deviated from the general tendency. We suppose the deviations were related to some local features of snow stratigraphy that attenuate fractionation despite shallow snow. In contrast to H2, d-excess showed significant correlations with SWE as well (Table 5). However, the strongest relationship was observed in 2019, when the SWE values largely corresponded to the distribution of snow depth. This

similarity was largely responsible for the significance of this relationship. In other sampling years, the linear influence of SWE on d-excess was not evident, despite the significance of the relationship (with the presence of outliers and heteroscedasticity).

**Figure 7.** Relationship of the d-excess and snow depth, SWE over the sampling years and basin parts (NP, CP, SP).

**Table 5.** Adjusted R2, slopes and intercepts for computed regression equations describing d-excess vs. snow depth and SWE relationships (ns = predictor variable not selected).


\* *p* < 0.1; \*\* *p* < 0.05; \*\*\* *p* < 0.01.

#### **4. Discussion**

The snow isotopic composition showed differences between open and forested areas, mainly in reducing d-excess and the slope of δ2H vs. δ18O regressions. The main predictor connected with snow enrichment in heavier isotopes was snow depth. The effect was especially evident in winter seasons with shallow snowpack.

The slope and intercept of δ2H vs. δ18O regressions corresponded to the values typical for the continental climate [41]. Additionally, the slope and intercept values agreed well with data previously obtained for the same region [31]. Negative intercept values and their high variability are typical for regions with a cold continental climate due to significant differences in the mechanisms of precipitation formation [41]. Under these conditions, the consistently low slope values and low d-excess values in the forest part most likely indicate the influence of sublimation.

The effect of canopy interception on the isotopic composition was observed not in the direct enrichment in heavier isotopes but through a decrease in d-excess. This pattern slightly differs from the existing studies in which the enrichment in heavier isotopes was correlated with increasing canopy density [26,27]. However, a recent study has shown that canopy density may not have a statistically significant influence on the spatial variability of δ18O [28]. Forest d-excess may be lower than in open areas (in [27], d-excess was not given but could be calculated from the average δ2H and δ18O values) or almost independent of canopy density [26]. We suppose that the lack of direct enrichment occurred due to the lower interception capacity of pine forests and low energy inputs for sublimation. Compared to fir and spruce forests, snow interception, and sublimation losses in pine forests may be up to 10% lower [42], despite the slightly longer deposition of snow on the canopy [43]. The work [27] noted slightly smaller interception effects on isotope ratios at the high-elevation transect, which the authors attributed to lower energy inputs. In the cold continental climate of southern Siberia, this factor may also limit sublimation.

The influence of snow depth on the isotopic composition was the most critical factor of both spatial and temporal variability. Previously, it has been shown that the upper and lower layers of snow are sensible to the changes in the isotopic composition due to the influence of melting, sublimation, and moisture exchange with the soil [6,20,28,44,45]. Vapor flux from the lower snow layers to the upper ones (due to the thermal gradient) usually does not significantly change the isotopic composition since all changes occur within the snowpack. These processes contribute to the homogenization of the isotopic composition between the individual layers during the winter [31,46]. In shallow snow conditions, bulk snowpack isotopic composition is much more sensitive since the upper and lower layers constitute a significant part of the snowpack. In the areas with increased snow accumulation (high snow depth), the fraction of "stable" snow within the snowpack is higher. Such snowpacks tend to be isotopically lighter. In moderate snow winter seasons, the isotopic composition-depth relationship weakens, but it is still recognizable through the d-excess decrease. In high snow winters, the correlation with snow depth almost disappears because a large amount of snow accumulates in the entire basin. The relationship between the isotopic composition and snow depth was earlier observed in alpine conditions [21]. The study [7] previously expected that samples with lower SWE should have been more strongly affected by melt or sublimation and, as a result, become isotopically heavier. However, they did not find a relationship between isotopic composition and SWE in the Canadian prairies. In our work, SWE also did not correlate with the isotopic composition. We argue that snow depth is a more significant factor since both higher density (due to wind exposure or sublimation) and higher snow depth (in a sheltered location) can produce the same SWE value. In terms of isotopic fractionation, these processes may result in different δ18O and δ2H ratios, which explains the lack of correlations.

The obtained results suggest that considering spatial variability and inter-annual differences in snow depth is necessary when planning observation strategies involving snow sampling for isotopic composition analysis. Additionally, spatial coverage is important because snowpack isotopic signatures exhibit a high spatial variability over a small scale (<100 m), limiting the usefulness of point samples to estimate an average isotopic composition of snow over a large area [7]. However, our findings also have limitations since we only estimated bulk snowpack isotopic composition, which implies some uncertainty in its evolution. To better understand the mechanisms of isotopic signal transformation, a joint analysis of the isotopic composition of snowpack, primary snowfall, and throughfall during the winter is preferable in further studies.

#### **5. Conclusions**

In the West Siberian forest steppe, the expected direct enrichment of forest snowpack in heavier isotopes was not observed. However, we found lower d-excess values and δ2H vs. δ18O regressions slopes less than those of GMWL in forested areas compared to open, which may indicate the influence of sublimation. Additionally, isotope ratios between open and forested areas maintained in both shallow and high snow winters. We found out that snow depth is the most critical spatial factor influencing the isotopic composition of snow. As snow depth increased, the δ18O and δ2H values of the entire snowpack decreased due to conservation within the snowpack and less influence of sublimation and moisture exchange with the soil. However, this pattern was only evident in winter seasons with below-average snow depth. In above-average snow winter seasons, significant amounts of snow accumulated at most of the sampling sites, which smoothed out the isotopic differences and contributed to the preservation of the δ18O and δ2H ratios within the snowpack.

**Author Contributions:** Conceptualization, D.P., N.M. and D.C.; methodology, D.P.; software, D.P.; validation, D.P., N.M. and D.C.; formal analysis, D.P. and N.M.; investigation, D.P., D.C., R.B., D.Z. and L.L.; resources, N.M. and D.C.; data curation, D.P. and R.B.; writing—original draft preparation, D.P. and N.M.; writing—review and editing, all authors; visualization, D.P. and R.B.; supervision, D.C.; project administration, D.P. and N.M.; funding acquisition, D.C. and N.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** Data collection was carried out within the State Assignment of IWEP SB RAS FUFZ-2021-0007. Analysis of the topography and land cover influence on the snow isotopic composition was funded by RFBR, project number 19-35-60006. Approbation of the methods for assessing the relationship between atmospheric parameters and landscape heterogeneity was supported by the Russian Science Foundation under grant 21-17-00135 (https://rscf.ru/en/project/21-17-00135/).

**Data Availability Statement:** The data and R code presented in this study are openly available in Mendeley Data. DOI: 10.17632/2g24d7pgyv.1.

**Acknowledgments:** The authors would like to thank Tatyana Papina, head of the Chemical Analytical Center at IWEP SB RAS, and researcher Alla Eirikh for performing the laboratory analysis of the snow samples.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
