Tradeoffs between Stand Volume and Understory Vegetation Diversity in Quercus wutaishanica Forests under Climate Change

Abstract

Natural forests play a crucial role in providing various ecosystem services, including timber production and biodiversity conservation. However, climate change and anthropogenic factors pose a severe threat to competing forest ecosystem services functions. Therefore, to optimize and sustainably utilize competing forest services, tradeoffs are often necessary. This study was conducted in Northwest China to explore tradeoffs aimed at improving the quality of Quercus wutaishanica Mayr natural forests under climate change conditions, focusing on stand volume, timber production, and understory vegetation diversity conservation. Data from 77 field surveys were used to construct a coupled model for stand growth, stand structure, and site conditions. Changes in understory vegetation species number (UVSN) with crown cover were quantified. These models and relationships can be used as tools to estimate tradeoffs. As stand density increased, single-tree volume decreased, whereas timber volume increased. UVSN increased and then decreased with increasing crown cover and was able to maintain a relative maximum at 0.5–0.65. Under the current climatic conditions, the optimum stand densities corresponding to 30, 40, 50, and 60 years were 1390, 1153, 1042, and 871 trees/ha, respectively, to maintain a high UVSN and adequate stand volume. When mean annual temperature rose, stand densities could be reduced to maintain high-quality timber. Although only two major services were considered, the tradeoffs presented in this study can inform future research to improve the quality of natural forests.

1. Introduction

Natural forests play a crucial role in ecology, climate, water resources, soil conservation, resource supply, and human well-being [1,2,3]. However, over-conservation measures, such as total logging bans and restrictions on human disturbance aimed at maximizing stand volume, have led to reduced efficiency in the provision of various forest ecosystem services (FESs) [4,5]; this has made it challenging to meet the increasing and diverse demands for services from society [6]. Consequently, traditional forest management approaches have been gradually revised to implement FES-oriented methods to enhance the quality and efficiency of natural forests [7,8].

Understory vegetation is an integral component of FESs and a vital indicator of forest health and stability [9,10]. It is strongly associated with numerous ecosystem service functions, including erosion control, carbon sequestration, hydrological regulation [11,12,13], natural tree regeneration [14], and biodiversity conservation [15].

However, natural forests still face numerous challenges, particularly under excessive protection measures; this includes high forest density leading to low understory vegetation diversity and cover [16,17,18], difficulties in seed regeneration [19], low timber quality, and poor resistance to snowstorms [20,21]. Therefore, it is crucial to establish multilayered forests that can provide multiple FES [22,23]. This includes a tradeoff between stand volume and understory vegetation species number (UVSN) with a certain level of timber quality [24,25]. However, whether there is a tradeoff between them and how this can be achieved remains controversial. Therefore, understanding and quantifying the responses of UVSN to tree and canopy densities is essential for guiding natural forest management [26].

Forest management is challenging due to the complex topography of natural forests located in mountainous areas [27]. Therefore, the development of models based on field measurements is particularly important to balance the competing services of stand volume, timber production, and understory vegetation diversity. In the context of a progressively warming climate, the introduction of climatic factors, including temperature and precipitation, into quantitative models can simulate the tradeoffs between changing FESs under multiple scenarios [28,29,30]. Stand and canopy densities directly affect forest UVSN, stand volume, and timber production functions [31,32,33]. Therefore, it is becoming increasingly important to balance these relationships by determining the optimal stand density in response to climate change [26].

Quercus wutaishanica is a widely distributed tree species in temperate deciduous forests in northern China. It is also a dominant species in the Liupan Mountain range. Large-scale logging was conducted in the 1950s. However, logging stopped in the 1980s with the establishment of nature reserves. Almost all the remaining forests are secondary forests formed by sprouting after logging. The main problems faced by forests are high stand density, low timber quality, low seed regeneration, and low understory plant species diversity. Therefore, there is an urgent need to improve the functions of multiple ecosystem services in secondary Quercus wutaishanica forests. To date, research on the ecosystem services of natural Quercus wutaishanica forests has mainly focused on carbon sequestration and hydrological regulation [34,35,36]. However, there has been relatively limited research on the response of UVSN to stand and canopy densities and integrated management guided by optimal stand and canopy densities to achieve multiple FESs [37,38].

This study aims to establish a coupled model of stand volume, single-tree volume, and UVSN responses to factors such as stand age and density based on field survey data. It also simulates the changes in these three functions to achieve their tradeoffs under climate change. This research can provide a theoretical basis for the comprehensive management of natural forests and how to efficiently utilize forest ecosystem services under climate change.

2. Materials and Methods

2.1. Study Area

The study area is located in the Liupanshan Nature Reserve (106°11′–106°32′ E, 35°20′–35°47′ N) in northwestern China (Figure 1). This region serves as a representative water conservation forest area in the northwest and plays a crucial role in soil and water conservation, and in ecological development in the Loess Plateau. The region is in a transitional zone from a warm temperate semi-humid zone to a semi-arid zone. The altitude ranges from 1700 to 2927 m. The mean annual temperature (MAT) decreases gradually from 6.5 ± 0.7 °C in the north to 5.8 ± 0.6 °C in the south (Figure 1). Mean annual precipitation (MAP) increased from 480 mm in the north to 548 mm in the south (Figure 1). The dominant soil type in the area is shallow and coarsely textured calcareous chestnut soil.

The forest cover of the Liupan Mountains region gradually increases from north to south, and the total cover of all forests is 65.4%. The main natural tree species are Quercus wutaishanica, Betula platyphylla, Populus davidiana, and Pinus armandii. Quercus wutaishanica forests account for 31% of the natural forest area. In the 1950s, the region experienced large-scale deforestation. Almost all the remaining natural forests are secondary forests formed by sprouting regeneration. The associated tree species primarily include Betula platyphylla and Populus davidiana. The understory shrubs consists of Viburnum schensianum, Cotoneaster multiflora, Spiraea blumei, and Crataegus pinnatifida. Meanwhile, the understory herbaceous plants include Epimedium brevicornu, Carex hancockiana, and Phlomis umbrosa.

2.2. Sample Plot Setting

Based on the main distribution range of Quercus wutaishanica in the study area and the characteristic of gradually increasing precipitation from north to south, sampling points were established for every 20 mm increase in precipitation. At each sampling point, four slope aspects were considered (defining the north direction as 0°, increasing clockwise, with azimuth angles 330°~30° as shady slopes, 30°~90°, 270°~330° as semi-shady slopes, 90°~150° and 210°~270° as semi-sunny slopes, and 150°~210° as sunny slopes) and three slope positions (upper slope, middle slope, lower slope). From 2021 to 2022, a total of 77 plots with an area of 20 m × 20 m or 30 m × 30 m were established and surveyed (Figure 1). In each tree layer survey plot, three shrub layer survey plots of 5 m × 5 m were established at the beginning, middle, and end positions along one diagonal line. Additionally, five herb layer survey plots of 1 m × 1 m were established at the four corners and the middle of the diagonal line of each tree layer survey plot.

2.3. Data Collection and Calculation

The elevation of the plot centers was measured using GPS. Slope, position, and aspect were determined using a compass. Crown cover was estimated as the ratio of the projected area of the canopy to the area of the sample plot. Tree cores were collected from all the trees in the plot using an increment borer, and the age of each tree was determined using the LINTAB 6 measurement system. The average age of the main canopy trees in the plot is used as the stand age for the plot. The tree density (trees/ha) was determined by dividing the total number of trees by the plot area. The height (H) and diameter at breast height (DBH) of all trees with a diameter ≥ 5 cm were measured. Mean H and DBH values for each plot were calculated using Equations (1) and (2), respectively.

$$H={\displaystyle \frac{1}{\mathrm{n}}}{\sum }_{\mathrm{i}-1}^{\mathrm{n}}{H}_i$$

(1)

$$DBH=\sqrt{{\displaystyle \frac{1}{\mathrm{n}}}{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{DBH}_i^2}$$

(2)

where $\mathrm{n}$ is the number of trees, and $H_i$ and ${DBH}_i$ are the H (m) and DBH (cm) of the ith tree, respectively.

The single-tree volume (${STV}_{trees}$, m³) of Quercus wutaishanica was calculated using the empirical equation developed by Meng et al. based on tree height and DBH (Equation (3)) [39]. The mean individual tree volume and stand volume (${SV}_{stand}$, m³/ha) of the forest were calculated by summing the individual tree volumes within each sample plot (Equation (4)).

$${STV}_{trees}=0.000057468552\times D_i^{1.915590}\times H_i^{0.9265972}$$

(3)

$${SV}_{stand}={\displaystyle \frac{\sum _{i=1}^n{V}_{{tree}_i}}{S}}$$

(4)

where ${STV}_{trees}$ is the single-tree volume, ${SV}_{stand}$ is the stand volume, $H_i$ is the tree height (m), $D_i$ is the tree diameter at breast height (cm), S is the plot area (m²), and n is the number of Quercus wutaishanica in the plot.

Based on the latitude, longitude, and altitude information of each plot, combined with the average precipitation and temperature data from six meteorological stations located in Guyuan City, Haiyuan County, Xiji County, Jingyuan County, Lund County, and Liupanshan Station uniformly distributed in the study area (data source: http://data.cma.cn/ (accessed on 3 October 2021)) from 1991 to 2019, and interpolated through Kriging interpolation using ArcGIS (Version 10.7, Redlands, CA, USA, Esri Inc.), the multi-year average precipitation and temperature values for each plot were obtained. The 77 sample plots investigated had a mean air temperature of 6.3 °C, a mean precipitation of 519 mm, and a mean slope orientation of 64°. In each of the 77 tree plots, we measured the plant species and numbers in three shrub plots and five herbaceous plots. To obtain the $UVSN$ for the tree plots, we calculated the mean values for the three shrub plots and five herbaceous plots and summed them. The specific calculation formula is as follows:

$$UVSN=\left({\displaystyle \frac{\sum _{i=1}^3{S}_i}{3}}\right)+({\displaystyle \frac{\sum _{i=1}^5{H}_i}{5}})$$

(5)

where $S_i$ represents the values from the shrub plots and $H_i$ represents the values from the herbaceous plots.

2.4. Construction of a Coupled Model of Stand Structure

By plotting a scatter plot of the dependent variable against the independent variable and selecting the highest points corresponding to the changes in the independent variable, these points of the dependent variable can be considered as being maximally influenced only by the independent variable. The line fitted to these selected points is the upper boundary line. Therefore, the upper boundary line can clearly determine the response function relationship of the dependent variable to the independent variable [40]. The modeling process was as follows, taking DBH as an example: (1) The correlation between DBH and each factor was determined. (2) The relative DBH was obtained by dividing the maximum mean DBH of all plots by the mean DBH of all plots. (3) A scatter plot of DBH versus the first factor, that is, the factor with the largest correlation, was constructed, and the upper boundary line was used to determine the response function $f\left(x_1\right)$ with the first factor. (4) To eliminate the effect of the first factor, all the data were divided by the value corresponding to the upper boundary line established by the DBH with the first factor. (5) The influences of the previous factors were eliminated individually until the influence function $f\left(x_n\right)$ of the last factor was determined. (6) The functions were multiplied and the model parameters were calibrated based on the measured data (i.e., DBH, height, corresponding stand density, annual precipitation, annual temperature, and aspect data) using 1stOpt (Version 1.5, Beijing, China, developed by Seven Dimensions High-Tech Co., Ltd.) to derive a coupled model of forest growth ($FS$) influenced by multiple factors. The functional form of the model is given in Equation (6).

$$FS=f\left(x_1\right)\times f\left(x_2\right)\dots \dots \times f\left(x_n\right)$$

(6)

where $FS$ is the mean H, mean DBH, or crown cover, $x_1, x_2,$...., $x_n$ are factors affecting H (m), DBH (cm), or crown cover, such as forest age and stand density, and n is the number of influencing factors.

Based on the mean H and DBH models and the expression of the binary volume equation, the mean single-tree timber volume (STV, m³) and stand volume (SV, m³/ha) models with the effect of the stand factor were constructed in the following model form:

$$STV=a\times D_p^b\times H_p^c$$

(7)

$$SV={D_p}^{b_0}\times {H_p}^{c_0}\times \left(d_0+N\times e_0+N\times {f_0}^2\right)$$

(8)

where $H_p$ and $D_p$ are the mean H (m) and DBH (cm). Their values were calculated using an established coupled model (Equation (5)) of stand mean H and DBH were affected by stand age, stand density (N, trees/ha), slope direction (SDP), MAT (°C), and MAP (mm). a, b, $b_0$, $c_0$, $d_0$, $e_0$, and $f_0$ obtained by re-rating using data calculated using the established coupled model (Equation (6)).

The cumulative breast height basal area (G, m³/ha) of the stand is an indicator reflecting radial growth and N and is correlated with crown cover. We established a relationship between G and the canopy by incorporating the developed DBH-coupled model to simulate crown cover variations under different scenarios. Therefore, by using the constructed DBH model in conjunction with N, we recalculated the G for each plot using Equation (9):

$$G={\displaystyle \frac{\pi {D_p}^2}{4}}\times N$$

(9)

2.5. Simulation of UVSN

Crown cover is a stand structure indicator that most directly affects understory vegetation species diversity. We constructed a response function relationship between crown cover and G. We then constructed a response function relationship between UVSN and crown cover. We embedded the constructed response function of UVSN to depression into the response function of crown cover to G and modeled changes in UVSN based on changes in G.

2.6. Principle of Tradeoffs

The tradeoff principle of this study is to achieve sufficient stand volume while appropriately adjusting the N to increase the potential UVSN. Meanwhile, this can take into account a certain amount of timber yield. The specific steps are as follows: (1) Based on the constructed coupled model and relationship of stand volume, individual tree volume, and UVSN, the changes in these three parameters with N under different forest age conditions were simulated. (2) The N was obtained corresponded to 90% of the maximum value of the changes in stand volume and UVSN. The N at each stand age served as the N interval for the upcoming tradeoff for each function. (3) The N intervals for each function overlapped, and the overlapped N intervals were at suitable density intervals that accounted for each function. (4) Given that the relationship between the individual tree volume and stand volume with N was negative, to take into account a certain timber yield, we selected the minimum density of the overlapping part of the stand volume and UVSN as the optimal density.

2.7. Model Evaluation

The model evaluation uses the coefficient of determination (R²) and the root mean square error (RMSE), which were calculated using Origin (Version 2021, OriginLab Corporation, Northampton, MA, USA). R² evaluates the fitting effect of the model on the modeling samples, while RMSE assesses the prediction effect of the model on the validation samples. Additionally, the Akaike Information Criterion (AIC) and the Bayesian Information Criterion (BIC) were calculated using R to assess model quality, considering both the goodness of fit and model complexity.

$$R^2=1-{\displaystyle \frac{\sum _{i=1}^m{(y_i-\widehat{y_i})}^2}{\sum _{i=1}^n{(y_i-\overline{y_i})}^2}}$$

(10)

$$RMSE=\sqrt{{\displaystyle \frac{\sum _{i=1}^m{(y_i-\widehat{y_i})}^2}{n-1}}}$$

(11)

where $y_i$ and $\widehat{y_i}$ represent the observed and predicted values of the forest structure index for the i-th sample, respectively, and $\overline{y_i}$ is the mean of all observed values. m is the number of samples used for parameter calibration or validation, and n is the number of samples used for model calibration or validation.

2.8. Scenario Formulation

Two scenario simulations were conducted to strike a balance between timber production and understory biodiversity conservation. In Scenario 1, we simulated the current climatic conditions, aligning precipitation, temperature, and slope with the average values observed in the study area. For Scenario 2, we introduced a specific change by elevating the temperature by 1.5 °C, in accordance with predictions outlined in the IPCC Fifth Assessment Report.

3. Results

3.1. Variation in Maximum Stand Density with Age

As trees age, stand density gradually decreases (Figure 2). The upper boundary line was used to estimate the maximum N at a given age. The maximum stand densities at the ages of 30, 40, 50, and 60 years were 2385, 2116, 1777, and 1368 trees/ha, respectively.

3.2. Response and Modeling of Stand Growth to Stand Structure

3.2.1. Response of Tree Height and DBH to Stand Structure and Site Factors

Forest stand growth is influenced by various factors (Figure 3). With the increase in tree age, mean H and DBH exhibited S-shaped growth patterns. During the 0–30 year period, their increase was gradual. In the 30–50 year period, the growth rate accelerated, and after 50 years, growth stabilized. They gradually decreased with an increase in N, experiencing a decline when N exceeded 1750 trees/ha.

Forest stand growth is influenced by age, density, and site conditions and climate (Figure 3). The mean H and DBH gradually decreased with an increase in slope direction, with an acceleration in decline when the slope exceeded 100°. They exhibited a parabolic relationship with mean temperature, rising when temperatures were below 6.4 °C and rapidly declining when temperatures exceeded 6.4 °C. They exhibited an increasing trend with increasing precipitation.

3.2.2. Construction of a Multifactor Coupled Model for Stand Growth Response

The upper boundary line functions of H and DBH in response to various factors were used to construct a multivariable coupled model to analyze the responses of average H and average DBH to various factors (

Table 1. Multifactor coupled modeling of stand growth.

Stand Structure	Multifactor Coupled Modeling of Stand Structure in Response to Various Environmental Factors	R2	RMSE	AIC	BIC	Number
H (m)	$H=33\times \left[\left(-1.47+{\displaystyle \frac{3.88}{1+{\left(\frac{Age}{2.25}\right)}^{0.58}}}\right)\times ((3.04\ast {10}^{-4})/(1+EXP((3.22\times {10}^{-4})\times (SD-914.93))))\times (779.66+3.25\times {MAP}^{4.24})\times (-10.73+13.27\times MAT-1.15\times {MAT}^2)\times ((-1.92\ast {10}^{-10})/(1+EXP(0.0035\times (SDP+80.49))))\right]$	0.72	1.36	373.8	410.7	(12)
DBH (cm)	$DBH=66\times \left[(0.006+(-0.003)/(1+{\left({\displaystyle \frac{Age}{41.64}}\right)}^{8.69}))\times (-35.31-890.14\ast {SD}^{0.46})\times (2.37\times {10}^5+0.29\times {MAP}^{2.36})\times (-0.0007-0.002\times MAT+0.0002\ast {MAT}^2)\times (0.08/(1+EXP(0.002\times (SDP+4174.47))))\right]$	0.80	3.77	358.0	394.9	(13)
Single-tree volume (m3)	$STV=2.24\times {10}^{-7}\times {DBH}^{5.06}\times H^{0.51}$	0.90	0.03	326.2	333.1	(14)
Stand volume (m3/ha)	$SV={DBH}^{1.53}\times H^{1.15}\times \left({-2.98\times 10}^{-3}+SD\times {1.48\times 10}^{-4}+SD\times {-2\times 10}^{-8}\right)$	0.85	25.45	338.5	347.2	(15)

). The complete dataset was divided into two parts. A total of 58 datasets were used for parameter calibration, and 19 datasets were used for model validation. The results demonstrated that the model performed well in terms of accuracy and precision (Figure 4).

The parameters of the formulas for individual tree and timber volumes were recalibrated using the constructed multifaceted coupled model for H and DBH (Table 1). The model demonstrated satisfactory performance (Figure 5).

The DBH calculated from the multifactor coupled modeling of DBH was introduced into Equation (8), and G was recalculated. The simulated values responded well to the measured values (R² = 0.73; Figure 6, Table 1).

3.3. Relationship between Understory Plant Diversity and Crown Cover

3.3.1. Changes in Crown Cover with Breast Height Basal Area

As G increased, crown cover ($CD$) exhibited an initial rapid increase (Figure 7). This was followed by a gradual slowdown when G reached 15–25 m³/ha. Beyond 25 m³/ha, crown cover ceased to increase. The relationship between G and crown cover was effectively captured by the function $CD=0.76(1-exp(-0.13G\left)\right)$, with a robust fit indicated by an R² value of 0.73.

3.3.2. Relationship between UVSN and Crown Cover

In the understory plant survey, we identified 188 plant species, including 72 shrubs and 106 herbs, belonging to 47 families and 98 genera. The most common shrubs were Cotoneaster multiflora Bunge., Spiraea blumei G.Don, Ostryopsis davidiana Decne., Smilax stans Maxim., and Viburnum betulifolium Batal. The common herbs were Carex hancockiana Maxim., Thalictrum aquilegiifolium var. sibiricum Linnaeus, Viola verecumda A. Gray, Sanguisorba officinalis L., and Discorea nipponica Makino.

Crown cover emerged as the dominant factor influencing UVSN compared with other structural indicators. As illustrated in Figure 8, UVSN exhibited a rapid increase with a progressive increase in crown cover, attaining elevated levels within the 0.5–0.65 range. Nevertheless, there was a rapid decrease in UVSN when the crown cover exceeded 0.65.

3.4. Tradeoff between Timber Production and UVSN

Compared to the maximum stand density achievable at each stand age under current climate conditions, achieving a tradeoff between stand volume and UVSN would require adjusting the maximum stand density at each stand age (2385, 2116, 1778, and 1369 trees/ha at 30, 40, 50, and 60 years, respectively) to 1390, 1153, 1042, and 871 trees/ha, respectively. Consequently, each growth indicator changes under the tradeoff scenario. Specifically, stand volume decreased by 7%, 10%, and 11% at 40, 50, and 60 years, respectively, while it increased by 2% at 30 years; UVSN increased by 4.79%, 10.16%, and 7.47%, except for a decrease of 10.81% at 30 years. Single-tree volume increased by 36.75%, 45.38%, 44.09%, and 37.43% at 30, 40, 50, and 60 years, respectively; this tradeoff demonstrates significant improvements in timber quality and UVSN (Figure 9).

After a 1 °C increase in temperature, compared to the tradeoff results under current climate conditions, stand volume and single-tree volume decreased by 33.33% and 44.05%, respectively, for all tree ages except UVSN. Additionally, compared to the maximum stand density achievable at each stand age after a 1 °C increase in temperature, achieving a tradeoff between stand volume and UVSN would require adjusting the maximum stand density at each stand age to 1666, 1158, 1000, and 866 trees/ha at 30, 40, 50, and 60 years, respectively. Similarly, under this tradeoff scenario, each growth indicator also changed. Specifically, stand volume decreased by 6%, 10%, and 15% at 40, 50, and 60 years, respectively, while it increased by 4% at 30 years; single-tree volume increased by 24.2%, 45.38%, 44.1%, and 45.45% at 30, 40, 50, and 60 years, respectively. UVSN increased by 8.58% and 8.06% at 50 and 60 years, while it decreased by −9% and −3% at 30 and 40 years, respectively (Figure 9).

Additionally, to restore the single-tree volume after warming to the level under the current climate condition tradeoff, the stand density needs to be reduced to 820, 700, 622, and 555 trees/ha at the ages of 30, 40, 50, and 60 years, respectively, under the tradeoff scenario after the temperature increase (Figure 9).

4. Discussion

4.1. Effects of Stand Structure and Climatic Conditions on Stand Growth

Stand growth is influenced by stand age and tree density [41,42,43]. Stand growth was described well using the logistic growth model in line with previous studies [44,45]. This finding aligns with the results of our study, showing that natural Quercus wutaishanica forests exhibit slow growth during the first 0–30 years, accelerated growth for 30–50 years, and slow growth beyond 50 years (Figure 3).

Increasing N exacerbates the competition for light, water, and nutrients, leading to reduced N [11,46,47]. In this study, stand growth began to decline rapidly when stand density exceeded 1750 trees/ha (Figure 3). This is consistent with the findings of Ahmad et al., where a power function effectively explains this result [40]. Although increasing stand density can enhance stand volume, it also adversely affects wood quality [48,49]. As N increased, single-tree volume decreased at an accelerated rate (Figure 3).

Stand growth is influenced by stand and climatic conditions [42,43,50,51]. Incorporating changes in temperature and precipitation into the model helped understand alterations in the structure–function relationship under climate change, enabling precise management [52,53]. Stand growth exhibited an initial increase, followed by a decline, as the mean annual temperature increased (Figure 3). Excessively low or high mean annual temperatures are detrimental to stand growth [54]. The optimal temperature for Quercus wutaishanica growth is approximately 6.4 °C. Considering the variations in the distribution of Quercus wutaishanica forests, precipitation conditions should also be incorporated into the model. As the mean annual precipitation increased, the growth of Quercus wutaishanica forests gradually increased (Figure 3). This relationship has been confirmed by numerous studies [42,43,55].

4.2. Effects of Stand Density and Stand Age on UVSN

UVSN decreases with increasing crown cover [24,56]. However, our study showed a different trend. UVSN initially accelerated with increasing crown cover, reaching a relatively high level when the crown cover was between 0.5 and 0.65 (Figure 8). Subsequently, UVSN gradually declined, corresponding to an N of 800–900 trees/ha. The reason for this phenomenon is that high stand density and crown cover favor the development of shade-tolerant plants while inhibiting the growth of shade-intolerant plants. As stand density increases, shade-tolerant plants gradually develop, resulting in an increase in the UVSN. However, when stand density and crown cover reach a certain threshold, shade-intolerant plants begin to die off. Therefore, with a further increase in stand density, the UVSN starts to decrease [57]. Therefore, UVSN reached its maximum value only when a balance was achieved between shade-intolerant and shade-tolerant species, in line with the findings reported by Zhang et al. [10]. Therefore, maintaining a crown cover within the range of 0.5–0.65 and an N of 800–900 trees/ha is the optimal approach to sustain a higher UVSN. However, when considering other dominant functions, such as timber production, which may vary with stand age and temperature, a fixed management approach may not be feasible [41].

4.3. Tradeoffs between Timber Production and UVSN

Among the various indicators of stand structure, N, stand age, and crown cover directly or indirectly influence timber production and the UVSN [58,59]. Stand density is the most decisive factor in forest management practices [41,60]. An excessively dense forest can result in increased stand volume but decreased single-tree volume. Meanwhile, excessive crown cover can lead to a decline in the UVSN [61]. In this study, we developed a generic tradeoff approach to maintain an optimal N corresponding to specific stand ages. Under current climatic conditions, the recommended stand densities are 1390, 1153, and 871 trees/ha for trees aged 30, 40, 50, and 60 years, respectively.

The model developed in this study incorporates climatic conditions—that is, temperature and precipitation—and slope direction as indicators, providing a technical means of balancing management based on different site types under climate change conditions. Compared with the model developed by Ahmad et al., which only considered stand density and stand age, our model included temperature, precipitation, and slope direction factors to make it more generalizable [62]. Therefore, as a more precise decision-making tool, the developed model can be used to determine the appropriate N at any stand age under climate change conditions, balancing the competing services of timber production and understory vegetation diversity conservation.

4.4. Application of the Model

MAP, MAT, and slope aspect are indicators that can be easily obtained from field surveys, facilitating the management of forests under climate change according to different site types. For a given site type, MAP and MAT can be interpolated from nearby meteorological stations or derived using remote sensing techniques, while stand age and slope aspect can be obtained through field surveys. Subsequently, using Equations (12) and (13) from this study, the simulated values of H and DBH can be calculated as they change with stand density. These results can then be incorporated into Equations (14) and (15) to derive changes in stand volume and single-tree volume with respect to stand density. Additionally, the relationship between crown cover and UVSN can be used to calculate changes in UVSN for that site type with varying stand densities. Finally, based on the tradeoff principles established in this study, the optimal stand density range that simultaneously satisfies timber production functions and the conservation of understory plant diversity for that site type can be determined.

4.5. Limitations of the Current Study and Directions for Further Research

The findings of this study have provided a reference for further research on the multifunctional management of forests with different tree species in different regions. However, the direct application of the model is limited because its age range is relatively low (15–66 years), which does not fully cover the rapid growth phase and the mature phase of stand development. Therefore, it is necessary to recalibrate the model parameters when applying it to older stands.

Soil drought is a significant factor limiting the growth of arid-land plants. Site factors such as elevation, slope, soil thickness, and soil texture have a significant regulatory effect on soil moisture. Therefore, these site factors and their impacts should be considered in future research.

Forests provide a variety of sustainable ecosystem services, including hydrological regulation, erosion control, seedling regeneration, and functions related to ecological quality. However, this study only considered timber and species diversity functions, and future research should incorporate more functions.

5. Conclusions

Stand growth is influenced by tree age, density, site conditions, and climate. H and DBH exhibit an S-shaped growth pattern with age, accelerating between 30 and 50 years and stabilizing after 50. When density exceeds 1750 trees per hectare, both H and DBH decline; they also decrease with increasing slope and temperatures above 6.4 °C, while increased precipitation promotes growth.

Based on the functional relationships of H and DBH in response to multiple single factors, a coupled model was developed to simultaneously assess the responses of DBH and H to stand density, tree age, mean annual temperature, mean annual precipitation, and slope. This model facilitates the tradeoff between timber production and the conservation of understory plant diversity under climate change.

Under current climate conditions, the tradeoff stand densities at stand ages of 30, 40, 50, and 60 years are 1390, 1153, 1042, and 871 trees/ha, respectively, allowing for a balance between timber production and understory plant diversity conservation. After a 1 °C increase in temperature, the tradeoff stand densities at these ages are 1666, 1158, 1000, and 866 trees/ha, respectively, necessary to achieve the same balance. Using the established model and functional relationships, more appropriate tradeoff decisions can be made to optimize both functions under varying tree ages and climate changes. However, when applying this model to other regions or tree species, recalibration of the parameters should be considered, and additional site conditions and forest service functions should be incorporated into the model.