Modelling the response of larch growth to age, density, and elevation and the implications for multifunctional management in northwest China

Plantations of Rupprecht’s larch (Larix principis-rupprechtii) have been widely established in the drylands of northwest and north China under traditional fast-growing plantation management strategies. These strategies and the long-term logging ban have led to over-populated stands with lower structural and functional stability, less economic benefit and higher water consumption. To guide the sustainable management of larch plantations, field surveys and historical data compilation were undertaken in the Liupan Mountains of northwest China. The main influencing factors (stand structure and site condition) and their effects on mean tree height, mean DBH and timber volumes were determined based on up-boundary line analysis. Tree growth models coupling the effects of tree age, stand density, and elevation were established. Both height and DBH markedly increased initially and then slowly with tree age, decreased with stand density, and showed unimodal change with elevation. The coupled growth models accounted for 72–78% of the variations in tree height, DBH and timber growth. Recommendations for future plantation management are: (1) prolong the rotation to at least 60 years to produce large-diameter, high-quality timber and maintain greater carbon stocks; (2) zone the target functions of stands by elevation; and, (3) reduce stand density for balanced supply of multiple ecosystem services. The growth models developed can predict growth response of larch plantations to density alteration under given ages and elevations, and assist the transformation from traditional management for maximum timber production to site-specific and multifunctional management with longer rotations and moderate tree density.


Introduction
Many ecological forestry programs have been implemented utilizing the important fast-growing native species, Larix principis-rupprechtii Mayr, in northwest and north China, where water shortages are problematic due to frequent droughts and severe soil erosion under sparse vegetative cover. These programs include the "Three-North Shelter Belt" begun in 1978, and the "Grain for Green" and "Natural Forest Protection Program" implemented since 2000 Wang et al. 2018), aimed at increasing forest cover for timber production and/or soil erosion control. However, many densely populated stands were created following traditional timber plantation management strategies which involved high initial planting densities to counteract low survival rates on difficult sites (Amazonas et al. 2018), minimal or non-existent thinning to maintain higher timber volumes (Sun et al. 2005), or maintaining dense vegetation cover to reduce soil erosion (Eric and François 1997). Furthermore, the traditional practice of the final harvest based on maximizing mean annual growth rates (Wang 2003;Duan et al. 2017) resulted a short 26-year rotation for larch plantations (Department of Forest Resources Management 2008), producing small to medium size timber of low unit price. Such management strategies, combined with the long-term logging ban policy implemented since 2000, have brought about many negative impacts, including low quality small timber of less economic value , significant water yield reduction due to high evapotranspiration (Yu et al. 2009), decreased structural stability against snow and storm damage (Hao et al. 2012), and restricted undergrowth and natural regeneration (Ahmad et al. 2018(Ahmad et al. , 2021. To solve these problems and address increasing societal demands for a balanced supply of various ecosystem services , China has issued a new national forest management plan for 2016-2050 (State Forestry Administration 2018) with the guiding principle of multifunctionality. While trade-offs among the conflicting demands on forest functions/services are unavoidable, this new plan is vital for achieving multifunctional forest management. Regardless of the trade-offs between forest cover and water management at a catchment level , or between commodity goods and ecosystem services at a stand scale (Sophie et al. 2017;Ahmad et al. 2018;Chen et al. 2018;Strengbom et al. 2018), timber production is always the most basic issue in trade-offs of various ecosystem services. However, these trade-offs consider every ecosystem service with equal weight over an extended area, neglecting the differential importance of various ecosystem services under diverse conditions of site, vegetation, and societal demand. Moreover, delivering ecosystem services of timber production involves individual trees with different economic values. There is a trade-off between the production of larger trees with higher unit price through lower tree density and the production of smaller trees with lower unit price through higher density. In fact, the production of high quality, large dimensional timber should be focused on fertile sites. Therefore, trade-offs between timber production and other ecosystem services should be regional-, site-, stand-and species-specific. A tree growth model reflecting the influence of multiple factors (e.g., site conditions and stand structure) is required to guide site-specific forest management. Most current models are empirical, e.g., the Richards equation (Richards 1959) and logistic equation (Kilian et al. 2014) which illustrate the effect of tree age, and the multiple linear regression model (Campoe et al. 2016) which considers the effects of meteorological factors. These empirical models are set out in a variety of mathematical forms but without a robust theoretical basis that affords the derivation of models with wide applicability (Valentine et al. 2005). In contrast, processbased models incorporate non-linear biological responses in a documented manner and provide rational explanations for causes of any growth changes (Coops et al. 2011). However, these models require numerous parameters to clarify growth mechanisms, which limits their fitting and application. Therefore, the development of hybrid models recognizing certain process-based variables are needed to guide timber production.
The establishment of fully or partially process-based models should be on the understanding of the dynamic effects of varying stand structure and the growth response to site variations. The interplay among stand structure, age, and density (Smith et al. 2001;Weiner et al. 2001) must be incorporated into the models as they are related to forest management. However, environmental factors and potential ecosystem services vary along gradients of site conditions, such as elevation affecting precipitation and temperature (Shi 2006), slope aspect and gradient influencing solar radiation (Méndez-Toribio et al. 2016;Renner et al. 2016), and soil depth determining the supply and storage capacities of water and nutrients (Tromp- Van et al. 2006). The integrated effects of these influencing factors on stand structure and site conditions must be incorporated into the models.
The objectives of this study were: (1) to understand and quantify the responses of growth indices (height, DBH, volume of average trees, and timber volume of stands) in larch plantations to site conditions (e.g., elevation) and stand structure (e.g., age and density); (2) to develop growth models using the effects of these factors; and, (3) to make technical suggestions based on the simulation analysis using fitted models for transforming the traditional timber management strategy to a multifunctional one.

Study area
The study area is within the small watershed of Xiangshuihe (106°12′-106°16′ E, 35°27′-35°33′ N) in the southern part of the Liupan Mountains, Ningxia, northwest China. The elevation ranges from 1900 to 2942 m a.s.l. The site has a temperate, semi-humid climate with a mean annual temperature of 5.8 °C, mean annual sunshine of 2100-2400 h, and mean annual precipitation of 676 mm, of which 73.3% occurs from June to September. The bedrock is composed of sandy mudstone and calcareous shale and the main soil types are grey-cinnamon (Cambisols) and alpine meadow soil of 0.4-1.0 m thick and rich in fragment content. The vegetation types are secondary broadleaved deciduous forests and meadows at higher elevations. The forest cover is 82.4%, of which 24.4% is plantations mainly of Rupprecht's larch.

Data collection and processing
Data from 273 plots were used in this study (Table 1) The variables included location, site factors (elevation, slope aspect and gradient, soil thickness), and growth parameters (tree age, density, height, DBH, canopy diameter and density, understory species and coverage). Timber volumes (V 0 , m 3 ) of average trees were calculated by the equation V 0 = 0.000057 × DBH 1.77 × H 1.125 including tree height (H, m) and DBH (cm) (Zhang et al. 2015a).
A preliminary analysis of stepwise linear regression of height and DBH response to influencing factors (age, density, elevation, slope aspect, slope angle and soil thickness) was carried out to identify the main factors for model establishment. It showed that the stand mean tree height was significantly affected by tree age (Age) and elevation (Ele) (H = 0.456Age + 0.004Ele − 8.297, R 2 = 0.708), while the mean DBH was significantly affected by tree age (Age), density (Den), and elevation (Ele) (DBH = 0.471Age + 0.007Ele − 0 .001Den −13.072, R 2 = 0.714). Therefore, age, density, and elevation were selected as the main factors influencing tree growth.
In this study, 217 plots were used for model parameterization by the upper boundary line (UBL) method and for model calibration, while other 56 plots were used for model validation. The reason for dividing the plot data into two groups was the even distribution of data within the variations of age, density, and elevation (Table 2).

Upper boundary line analysis
To develop the integrated models including the effects of all main influencing factors (age, density, and elevation) on tree growth (height and DBH) based on field data, the response trend and corresponding optimal function type of the growth indices to each individual influencing factor must be determined. To achieve this, UBL analysis was performed to derive the optimal function type through minimizing the effects of other factors (Mimra et al. 2014;Irantzu et al. 2015;Matsushita et al. 2015).
For UBL analysis, the variation ranges of independent variables (i.e., the influencing factors) were first divided into segments. The UBL data of the dependent variables were then selected for each segment if their values were higher than the mean plus one standard deviation within each segment (Ren et al. 2016). However, the data were not uniformly distributed among the independent variables, causing more scattering or insufficient data in a few segments. In this case, some UBL data within the segments containing very few data were artificially selected, or some redundant data within one segment containing excessive data were removed.

Establishment of height and DBH models
The structure of tree growth model is assumed to be a continuous multiplication of the response functions to all main influencing factors described in Eq. 1. Such a coupling method has been widely used for developing process-based models (e.g., Fang et al. 1999;Waring et al. 2011), including tree growth models (Matsushita et al. 2015).
where TGI represents specific tree growth indices (mean tree height or DBH of each plot); TGI max is the specific theoretic maximum value of TGI; f(Age), f(Den), and f(Ele) designate growth response functions to the influencing factors of age, density, and elevation, respectively. These functions were designed to vary from 0 to 1. The concrete forms and initial parameter values of these functions were determined using the UBL method and a decoupling process. For model establishment, the values of the constant TGImax were derived from Chinese Flora (Zheng and Fu 1978) and other literature describing natural and mature larch forests at ages close to or at their TGI max . All TGI were divided by TGI max to create dataset 1 of relative TGI values for further analysis. A decoupling process was then employed to determine the function types of tree growth response to each influencing factor (i.e., f(age), f(den) and f(ele)): (a) A scatter diagram of the relative TGI values in dataset 1 plotted against tree age was used to derive the function type of f(age) through UBL analysis; (b) The relative TGI values in dataset 1 were divided by the corresponding values of each data point on the UBL curve of f(age) to get dataset 2 for decoupling the effect of tree age. Data in dataset 2 were then plotted against tree density to derive f(den) through the UBL analysis; (c) Using the same procedure, dataset 3 was obtained to decouple the effect of tree density for deriving the f(ele). Finally, the framework of the coupled TGI model was established by multiplying the TGI max , f(age), f(den) and f(ele). These processes formed the model parameterization which then needed to be further calibrated and validated using original inventory data from forest plots.

Calibration and validation of timber volume models
The timber volume of average stand trees (V tree , m 3 ) was calculated by V tree = 0.000057 × DBH 1.77 × H 1.125 , using the stand mean DBH (cm) and height (H, m), while the stand timber volume (V stand , m 3 ‧ha −1 ) was calculated according to V stand = ∑V tree . A model of the timber volume of average stand trees (V tree , m 3 ) based on the stand mean TGI was established by incorporating the tree height model and DBH model into Eq. 2: where the model parameters of a, b and c were fitted using original plot data. Similarly, the model of stand timber volume was established according to Eq. 3: where N represents the tree density (ind.·ha −1 ), and the model parameters of a 0 , b 0 , c 0 and d 0 were also fitted using original plot data.

Limitation of the maximum tree density
For growth simulation of the even-aged and one-layered plantations, the limit of theoretical maximum tree density (N max , ind.·ha −1 ) with changing tree size (Gadow and Hui 1999) must be considered. This means that the density of live trees cannot exceed N max during model simulation. Therefore, a model was derived using the UBL analysis of plot data to show the N max variation with the mean DBH of trees.

Mean annual timber volume growth rate
To determine the final harvest age in traditional timber management, the maximum mean annual growth rate (MTGR , m 3 ·ha −1 ·a −1 ) is computed as: where MTGR of larch plantations was calculated based on model simulation results.

Relative effect of individual influencing factors on timber volume
Although the effect of individual factors on timber volume can be reflected from the UBL, it is difficult to quantify and compare the effect of certain single factors based on the variation of UBL alone. Hence, the relative effect of each individual factor was calculated under given values of two other factors as below, allowing variation within the range of 0-1 for easier comparison (Li and Liu 2011;Kılkış 2016): where RETG i is the relative effect on timber volume of the indices of V tree , V stand and MTGR by the influencing factor i (age, density, elevation), TG i is the timber volume simulated by the models under the influence of factor i with given values of two other factors, and Max (TG i ) is the function for defining the maximum value of TG i . For an example of how to evaluate the influence of elevation to timber growth, the relative volumes of average stand trees and stands (RETG ele ) are divided into three categories by two thresholds (0.6 and 0.8) which are given artificially. When the RETG ele is 0-0.6, 0.6-0.8 and 0.8-1, the potential production of timber volume is defined as least suitable, suitable, and optimal, and the corresponding elevations can be defined as least suitable, suitable, and optimal for timber production.

Statistics
Statistical analyses used Microsoft Excel 2010; the fitting of UBL equations and models with 1stOpt 1.5 based on the nonlinear least squares method, and the figures drawn with Origin 8.0.

Effects of individual factors on height and DBH
The reported maximum height and DBH of larch trees are 30 m at 80 years and 100 cm (no age recorded, according to the limited literature) (Zheng and Fu 1978;Zhang et al. 2015a). However, these studies did not examine tree growth and influencing factors simultaneously, and there is little consensus among larch inventory results, and therefore difficult to make a scientific conclusion. The potential maximum tree height and DBH were preliminarily assumed to be 33 m and 110 cm, respectively, when all influencing factors were optimal, and these values were used in later analysis.
The effects of individual factors on height and DBH were obtained through the UBL analysis (Fig. 1). Stand mean height rapidly increases with age up to 30 years. It reaches 68% of the theoretical maximum at 50 years of age (Fig. 1a1). The growth pattern of stand mean DBH increases with age like that of height but with smaller increments after 30 years. Stand mean DBH at 50 years only reach 26% of the theoretical maximum (Fig. 1a2).
Both the stand mean height and DBH decrease with tree density, but at varying rates in different densities. Height decreases slowly with density at up to 1200 ind.·ha −1 but thereafter decreases more rapidly, with a 25% reduction compared to its maximum. DBH is more sensitive to density, decreasing at a declining rate within the range of ˂ 3600 ind.·ha −1 , leading to a final 60% reduction.
The effect of elevation on growth follows a unimodal curve for both stand mean height and DBH, although the curve for height is less distinct. Peak height and DBH are found at elevations of 2350 and 2400 m, respectively. At 1800 m, height and DBH are the smallest, 40% and 45% lower than their maximum, respectively.

Establishment of coupled tree growth models
The framework construction and parameterization of growth models were developed by coupling the response functions of UBL with individual factors identified in Fig. 1 as below: where Age, Den and Ele represent the variables of tree age (a), density (ind.‧ha −1 ), and elevation (m), respectively. The letters a to g are model parameters to be fitted using plot data of group 1 (217 plots). The fitted model parameters for height and DBH are given in Table 3, with an acceptable level of prediction accuracy (R 2 = 0.76 and P < 0.01 for height; R 2 = 0.78 and P < 0.01 for DBH). The growth models were validated using the plot data of group 2 (56 plots), showing a more acceptable level of prediction accuracy (R 2 = 0.80 and P < 0.01 for height; R 2 = 0.85 and P < 0.01 for DBH) (Fig. 2).
The timber volume models for both individual trees (average stand trees, Eq. 2) and stands (Eq. 3) were fitted using the first group data of 217 plots with an acceptable level of prediction accuracy (R 2 = 0.77 for individual trees; R 2 = 0.72 for stands) (Table 4) and validated using the second group data of 56 plots with a more acceptable level of prediction accuracy (R 2 = 0.92 and P < 0.01 for individual trees; R 2 = 0.86 and P < 0.01 for stand), indicating a high applicability of the timber volume growth models (Fig. 3).

Simulated variation of timber volume with main influencing factors
From the UBL in Fig. 4, the maximum density (N max ) can be high when mean DBH is ˂ 10 cm. When the DBH is > 10 cm, N max decreases quickly with DBH up to 20 cm and then gradually thereafter. When the mean DBH is 10, 20, 30 and 40 cm, the corresponding N max is 4218, 1146, 535 and 311 ind.·ha −1 , respectively.
Based on the fitted growth models, the response of timber volume for stands and for average stand trees to increasing density at different ages (10, 20, 30, 40 and 50 years) and elevations (1800, 2100, 2500 and 2800 m) is illustrated in Fig. 5.
Timber volumes of both stands and average stand trees increases with age up to 50 years, regardless of density and elevation (Fig. 5). Taking an example of a fixed density of 500 ind.‧ha −1 , the mean timber volume of average stand trees for elevations of 1800, 2100, 2500 and 2800 m at 10, 20, 30, 40, and 50 years of age are 0.0031, 0.044, 0.156, 0.314 and 0.47 m 3 , respectively, while the corresponding means of Table 3 Parameters fitted for height and DBH growth models of Eq. 6 The letters a to g are model parameters to be fitted using plot data of group 1 (217 plots) Timber volumes of average stand trees decreased almost linearly with tree density (Fig. 5). Since densities over 1000 ind.·ha −1 do not exist at 40 and 50 years of age, only data covering 10, 20 and 30 years of age and elevations of 1800, 2100, 2500, and 2800 m were used to calculate the volume means of average stand trees at 500, 1000, 1500 and 2000 ind.·ha −1 , and they are 0.072, 0.061, 0.049 and 0.027 m 3 , respectively. In contrast, the stand timber volume increases with rising density up to a certain level, because either the maximum volume or the maximum density is reached, and then decreases gradually. The means of stand timber volume at densities of 500, 1000, 1500 and 2000 ind.·ha −1 are 57.1, 73.1, 78.3 and 58.2 m 3 ·ha −1 , respectively.
Timber volumes of both stands and average stand trees show unimodal change with rising elevation. For example, for a density of 500 ind.·ha −1 , the volume means of average stand trees (10, 20, 30, 40 and 50 years) at elevations of 1800, 2100, 2500 and 2800 m were 0.06, 0.19, 0.33 and 0.21 m 3 , respectively. The corresponding means of stand timber volumes were 50.25, 122.5, 176.34 and 126.28 m 3 ‧ha −1 , respectively. Figure 6a shows the variation of relative mean annual growth rate of stand timber volume with age. According to the principle for determining the final harvest age in traditional plantation management, all trees should be harvested when mean annual growth rate of stand timber volume reaches its maximum at 50 years (Fig. 6a). The relative timber volumes for both stands and average stand trees increase with age and follow logistic curves, but they are kept at a high level at the age of 60 years (Fig. 6b).

Implications of simulated growth for forest management
The effects of elevation on timber volumes of both stands and average stand trees are unimodal. The optimal elevation is 2170-2730 m for stand timber volume and 2230-2700 m for average stand trees. Therefore, the common optimal elevation for both higher stand volume and large-sized single trees is 2230-2700 m, as indicated by the light grey covering in Fig. 6c. The least suitable elevations for timber volume of stands and average stand trees are below 2050 and 2100 m respectively. This means that the least suitable elevation is < 2050 m for timber production of both individual trees and stands, as indicated by the dark grey covering area in Fig. 6c. Therefore, the remaining elevations (2050-2230 and 2700-2800 m) can be defined as the suitable elevation for timber production.

Comparison of DBH distribution among tree density ranges
Data of plots located within the optimal elevation of 2230-2700 m and with tree ages 25-35 years were selected to evaluate the density effect on DBH distribution when DBH has intervals of 3 cm (Fig. 7). The range of DBH distribution for both tree number and timber volume become narrower with tree density. If defining the borders of DBH distribution as the position of zero, then DBH distribution is 6-33 cm for the index of tree number and timber volume at a density of 500-1000 ind.·ha −1 . The corresponding DBH ranges are 6-30, 6-27 and 9-21 cm at densities of 1000-1500, 1500-2000 and 2000-2500 ind.·ha −1 , respectively. The frequency of trees with a DBH > 30 cm (volume > 0.43 m 3 ) declines sharply with tree density, i.e., they are 3.9% (8.3%), 0.1% (2.1%), 0% (0%) and 0% (0%) within the four tree density ranges, respectively. The DBH distribution concentrates at the ranges of 18-24 cm (accounting for 71.1% of total DBH classes), 15-21 cm (75.8%), 15-21 cm (74.4%) and 12-18 cm (82.0%), with a corresponding mean DBH of 21.2, 17.6, 16.7 and 15.2 cm, respectively. The DBH class with the peak value for timber volume distribution is one DBH interval larger than that for tree number distribution.

Age effect on tree growth
Changes in the relative growth rate with age are greater for height than for DBH within the ages of 10-50 years in this  (Fig. 1). This is physiologically rational but contrasts with a study on larch plantations 30-80 years of age in north China (Zhang et al. 2015a). The absence of the fast-growing period under 30 years for height in the latter may explain the difference between the two studies. The simulated relative timber volume for stand is greater than that for average stand trees before 60 years (Fig. 6b), but opposite after 60 years. This may be explained by the increasing tree competition with age, leading to the death of some suppressed trees but a continuous growth of remaining trees, especially the dominant and sub-dominant ones (Smith and Long 2001). The volume loss due to the death of suppressed trees will lead to a reduced growth of stand timber volume.

Density effect on tree growth
Forest growth is associated with and affected by fundamental changes in stand structure, including canopy reorganization, size-class differentiation, and self-thinning (Smith and Long 2001). Among the numerous stand structure parameters, density is the easiest one to alter for regulating timber growth of both stands and single trees. DBH is more susceptible to density variation than height (Fig. 1), which conforms to the general principle. According to size-density relations, as shown in Fig. 1, stand mean DBH decreases with rising tree density (Yang et al. 2002), while stand mean height is stable within the medium density range, and only sensitive to density in overly dense or sparse stands (Meng 2006).
In this study, the decreasing trend of timber volume for average stand trees with rising density (Fig. 5) is similar to that of DBH (Fig. 1); while the stand timber volume varied in a unimodal pattern (Fig. 5) due to the varying counteracting effects of increasing tree numbers and decreasing timber volume of average stand trees. The density effect on volume of average stand trees is a result of the intra-specific competition for space and carbon allocation (Weiner et al. 2001;Binkley 2004;Matsushita et al. 2015). Extensive studies have indicated that the occupation of space is linked with above-ground competition for light among dominant and suppressed trees, especially on fertile sites, whereas underground biomass is determined by below-ground competition for nutrients and water among all trees in a stand, especially on infertile sites (He et al. 2000;Pretzsch et al. 2010).

The effect of elevation on tree growth
The variation in volume growth with elevation within 1800-2800 m in this study was unimodal for both average stand trees and stands (Fig. 6c). Several studies on the montane forests in the dry northwest China reported that tree growth firstly increases with elevation and then decreases (Zhang et al. 2015b;Yang et al. 2017;Wang et al. 2019). However, other studies from humid tropical regions were in contrast to this, indicating a continuous decrease of tree growth with rising elevation (Coomes et al. 2007;Moser et al. 2008;Malhi et al. 2010). This might be explained by the different dominant climate factors and their vertical distribution in mountainous areas which are represented by thermal and precipitation gradients (Irantzu et al. 2015). In southern Ecuador, for example, where abundant annual precipitation of 2000-2300 mm within elevations of 1000-2000 m, water was not a limiting factor, whereas tree growth was influenced by temperature. In contrast, in the semi-humid region of this study, the limiting factors to tree growth were water availability which increases with elevation due to increasing precipitation and with low temperatures which decrease with elevation.

Management strategies for larch plantations
The timber volume growth response to density differs between stand and average stand trees. There is a problem between the goals of maximizing timber production and maximizing large-sized timber of individual trees from an economic perspective. Traditional timber plantation management are characterized by the pursuit of maximum mean annual growth of timber volume and the neglect of large dimensional timber which brings higher economic profit. This unavoidably leads to a short rotation clear cut (Qi 1992;Ou et al. 2018), e.g., 26 years for larch plantations (Department of Forest Resources Management 2008), and a more frequent disturbance  to forest/ vegetation cover, forest soil, and ecosystem services. However, preventing soil erosion, ensuring quality water supply, and maintaining high biodiversity and carbon stocks are also important ecosystem services to sustain the development of the study region. Thus, the function of high-quality timber production should be viewed as the dominant ecosystem service at some fertile sites, but not on infertile sites where a dominant ecological service (e.g., water yields in a dryland region) should be pursued. From this viewpoint, the traditional timber plantation management should be updated to the site-specific and multifunctional management (Duncker et al. 2012).
The timber volume growth rate varies significantly with elevation for both stands and average stand trees (Fig. 6c). Thus, the elevation range (1800-2800 m) in this paper was zoned for site-specific and multifunctional management. The 2230-2700 m range is optimal for timber growth of both stands and average stand trees, especially for the production of large-size timber, i.e., more weight should be given to timber production than to other ecosystem functions with higher commercial values but with a relatively high level of ecological supply. On suitable sites within the elevations of 2050-2230 and 2700-2800 m, equal weight should be given to timber production and other ecosystem services, or certain services which have optimal potential at these elevation ranges can receive more weight in the process of trade-offs among various ecosystem services. On least suitable sites with a drier climate and less fertile soils at elevation < 2050 m, the dominant functions should be noncommercial ecosystem services such as erosion control and/ or regulating water yields, as the study region is a source of water important for the surrounding drylands (Tian et al. 2021).
Based on the timber growth with age observed in this study (Figs. 6a and 6b), the previously proposed standard rotation of 26 years which causes frequent harvest disturbance should be prolonged to 60 years or more for attaining higher timber volumes and less disturbance to the stand structure and ecological services (e.g., carbon sequestration and soil erosion control), and producing high-quality, large sized timber for higher economic income due to greater unit price of timber. In fact, the proper rotation length for an optimal economic return should be based on the trade-off between timber production and other ecosystem services, as demonstrated in a deciduous forest in Turkey (Hayati et al. 2015).
Besides rotation length, keeping an optimum stand density is also important to balance competing ecosystem services and the conflict between timber income and these services. There are two aspects to be considered in determining an optimum density with multifunctionality management. The first is the economic trade-off between the higher price of large-size timber and the higher total stand volume. For example, a relatively low density of 500-1000 ind.·ha −1 is preferred for producing large-size timber (see Fig. 7). The second aspect is the trade-off among competing ecosystem services. For example, a canopy density of 0.7 and a corresponding stand density of 1255 ind.·ha −1 at 30 years were suggested to balance competing objectives of timber production and understory diversity protection based on a study in the same area (Ahmad et al. 2018). When there is a disagreement between economic and ecological values, a compromise is required, e.g., the multifunctional stand density at a given tree age and site quality must be determined through several trade-offs among the benefits of ecological services and stand stability. Such a comprehensive decisionmaking procedure for multifunctional forest management was proposed by Tian (2019).
The above suggested forest management measures, i.e., longer rotation and optimum stand density, is based on the calculated timber growth models fitted in this study. However, in practice, the optimum rotation length and stand density are still dependent on the maximization of economic profit, which is related to the market timber price, tax, and subsidies (Van et al. 1995;Kaiser et al. 2013). The current timber price in China does not differ significantly according to timber size and quality, possibly due to the shortage of large dimensional wood and related processing. This situation must be changed with the development and implementation of appropriate policies and a well-developed wood products market. Comprehensive consideration of the optimal configuration of stand structure, the corresponding ecosystem services, and the relative economic benefit after a reformed payment for these services should be considered simultaneously. An accurate DBH distribution model predicting the long-term DBH class variations is also needed to drive the development of a timber dimension and ecosystem services targeted management plan.

Conclusion
From this larch plantation growth study and multifunctional management analysis undertaken in the semi-humid Liupan Mountains in northwest China following conclusions could be derived: (1) The main influencing factors on tree growth are age, density, and elevation in the study region. The responses of tree height and DBH growth to these three factors were understood and quantified. They present a positive, negative, and unimodal relation, respectively. (2) The timber volume growth models for stands and average stand trees were established through coupling the effects of main influencing factors, and they can well predict the timber volume response to changing stand structure (age, density) and site condition (elevation). (3) Suggestions were derived for improving current forest management: (a) The previously proposed standard harvest rotation of 26 years for larch timber plantation is too short, and should be prolonged to at least 60 years for fully utilizing the growth potential of largesize timber and for lowering the disturbance to multiple ecosystem services, such as carbon sequestration and soil erosion control; (b) The whole elevation range (1800-2800 m) studied can be divided into three multifunctional zones with different combinations of dominant and other ecosystem services, i.e., the fertile sites (within the elevation range of 2230-2700 m) mainly for producing high quality timber, the least suitable sites (< 2050 m) mainly for water yield and/or soil erosion control, and the remaining sites (within 2050-2230 and 2700-2800 m) for both timber production and the supply of other ecosystem services; (c) The over-populated stands should be thinned to an optimal density range to enhance the forest stability, economic productivity, and