The study area is located in the northernmost and highest latitude region of China (50°07′02”-53°33′42”N, 121°10′53”-127°01′21”E; Figure 1). The terrain is predominantly characterized by hilly and low mountainous landscapes, with an average elevation of approximately 573 m above sea level, rising to a maximum of about 1528m. The region experiences a typical continental monsoon climate, characterized by a mean annual temperature of -2.45°C, an average annual precipitation of approximately 450mm, and a frost-free period lasting 80–100 days. The dominant soil type is dark brown loam, accompanied by smaller proportions of white slurry and meadow soils. The dominant tree species in the region include Larix gmelinii, Pinus sylvestris, Betula platyphylla, and Populus davidiana.

The data used for this study were obtained from the fixed plots of the 7^th National Forest Resource Continuous Inventory (NFI) of China. The plots were established on a systematic 4-km grid, each square in shape and 0.0667 hectares in area. According to the Technical Provisions for National Forest Resource Continuous Inventory (State Forestry Administration, 2004), a total of 243 NFI plots were selected from the database, ensuring that L. gmelinii accounted for more than 65% of the total stand basal area. For each plot, the recorded variables included geographic coordinates, site factors (e.g., elevation, slope aspect, slope gradient), soil factors (e.g., soil type, soil thickness), stand factors (e.g., stand origin, canopy density), and tree factors (e.g., species, diameter) (Table 1). In accordance with the NFI protocol, all trees with a diameter at breast height (DBH) ≥ 5cm were measured and recorded. The measurements included tree species, DBH, vitality status, and relative positions (radius and azimuth).

The living biomass carbon of each tree, partitioned into root, stem, branch, and leaf components, was estimated using the allometric biomass models developed by Dong (2015) and the carbon concentration coefficients reported by Jia (2014). The adopted carbon fractions were 0.52 for roots, 0.47 for stems, 0.51 for branches, and 0.52 for leaves, respectively. The component-specific carbon stocks were first calculated at the plot level and subsequently converted to a per-hectare level. The corresponding climatic variables for each plot, including mean annual temperature (T) and mean annual precipitation (P), were extracted using the ClimateAP model. In addition, the De Martonne aridity index [MAI = P/(T + 10)] was employed to quantify the site-level aridity condition. Generally, higher MAI values correspond to lower degrees of aridity. Statistics analysis showed that the MAI across the study area ranged from 39.5 to 95.7 mm/°C, with a mean of 61.7 mm/°C and a standard deviation of 10.6 mm/°C (Table 1), reflecting a wide range of hydroclimatic conditions.

The parameters for Equation 3 were estimated using linear quantile regression (rq function in “quantreg” package) at a 0.99 quantile. In contrast, the parameters for the remaining equations were determined through weighted nonlinear least-squares regression (nls function in “stats” package). The weighting functions were defined as the reciprocal form of the corresponding model function (Zeng et al, 2021). Model performances were evaluated using three statistical criteria: adjR² (Equation 14), MAE% (Equation 15), and RMSE% (Equation 16), respectively. A larger adjR2 value, together with smaller MAE% and RMSE% values, indicates better model performances.

where yi and y^i are the i-th observed and estimated values; y¯ is the mean value; n is the total number of observations; k is the number of parameters.

For a specific stand, SCI can be predicted using Equation 2, which was assumed to remain constant throughout the simulation period. The temporal dynamics of key stand variables (e.g., TH, DBH, BAS, NUM, VOL, and CAR) can be predicted following the procedure below.

5. Estimate N2 (Equation 18) at time t₂ using the definition of BAS ( BAS2=N2(π/40000)Dg22). and the derived Dg2 from step 4:

6. Predict total VOL (or CAR) using the obtained BAS (Equation 8) and TH2 (Equation 9);

7. Estimate carbon stocks for different components (leaf, branch, stem and root) using Equation 9;

8. Allocate the total stand density (N) into different diameter classes ( Ni) based on Step 5 and Equation 10;

9. Calculate the volume for various assortments (large-, medium-, small-, short-, fuel-wood and bark) and the carbon stocks for different components (leaf, branch, stem and root) using the PWA model, wood density and the estimated Ni values from Step 8.

Here, N₁ and N₂ denote the stand densities at time t₁ and t₂, respectively; Dg₂, BAS₂, and SDI2 represent the mean DBH, basal area, and stand density index at time t₂.

The parameter estimates and goodness-of-fit statistics for the FGYM models were summarized in Table 3. For the SCI model, the adjR2, MAE%, and RMSE% values were 0.3347, 17.23%, and 21.03%, respectively. Based on Equation 2, the frequency distribution of SCI values at the plot level was illustrated in Supplementary Figure 2, showing a mean value of 13.91 m and a standard deviation of 2.93 m. The slope of the SDI model (Equation 3) was significantly shallower (-1.2338) than that of Reineke’s self-thinning law (-1.605; Reineke, 1933). At the plot level, the mean and standard deviation of SDI values were 587.60 ± 282.02 trees hm^-2 (Supplementary Figure 3). Although the SDI dynamic model (Equation 5) exhibited relatively lower accuracy than the other models, with an adjR2 of 0.1552, MAE% of 36.53, and RMSE% of 44.28, it still effectively captures the relationship between stand age and SDI, as illustrated in Figure 3.

The accuracy of the base BAS model (Equation 6) was relatively low, with an adjR2 of 0.0407. However, a significant improvement in adjR2 (0.8815) could be observed when both SCI and SDI were incorporated into the model (Equation 7). Further inclusion of climate-related variable (Equation 12) and stage-related variable (Equation 13) enhanced the model accuracy by approximately 1.03% and 2.54% in adjR2, respectively, when compared to the stand-related BAS model (Equation 7). Correspondingly, notable reductions in MAE% and RMSE% were also achieved for Equation 12, 13, showing decreases of 3.53% and 11.57% in MAE%, and 4.54% and 12.56% in RMSE%, respectively. As illustrated in Figure 3, substantial differences in the estimated BAS were evident among the four BAS models.

The VOL model exhibited robust statistical performance, with an adjR2 of 0.9806, MAE% of 6.00 and RMSE% of 7.60, collectively indicating a high level of model accuracy. Figure 4 presented the residual distributions and the correlations between the measured and predicted values for the TH, SDI, BAS, and VOL models. These graphical assessments confirmed that all models satisfied the fundamental prerequisites for biometric modeling, e.g., the assumption of homoscedasticity.

The developed FGYM enables the simulation of (i) stand structure attributes, including Dg, TH, BAS, NUM and SDI); (ii) volume yields across different timber assortments (large-, medium-, small-, short-, fuel-wood and bark); and (iii) carbon stocks in various tree components (leaf, branch, stem and root) and end-use categories (commercial- and noncommercial- wood) under different combinations among MAI, SCI and stand age (see details from Appendix). Comparative analyses between the Stage-FGYM and the Base-FGYM revealed significant divergences in both stand volume and carbon stock trajectories (Figure 5). For a representative stand (age=100 years, SCI = 16m), the Base-FGYM predicted baseline estimates of 126.42 m³ hm^-2 (volume) and 57.80 t hm^-2 (carbon), while the Stage-FGYM simulations indicated clearly stage-dependent variations: at Stage 2 (Dg=16 cm), volume and carbon stocks increased by approximately 2.96% and 3.11% respectively, whereas they decreased by 15.02% and 15.70% at Stage 3 (Dg=20 cm), demonstrating the model’s sensitivity to ontogenetic transitions and developmental phase shifts.

Table 4 provides an example illustrating the potential effects of different MAIs (50, 60, and 70), SCIs, and development stages on stand characteristics, volume, and carbon stock when stand age is fixed at 80 years. For a specific MAI, Dg, BAS, VOL, and CAR all exhibited significant increases with rising SCI, irrespective of development stages, whereas NUM was strongly correlated with the developmental stage. In contrast, for any specific SCI, Dg, BAS, VOL, and CAR showed marked declines as MAI increased. Across all combinations between MAI and SCI, predicted values of Dg, BAS, VOL, and CAR at Stage 1 were notably larger than those obtained at Stage 0 (baseline, without development stage information), while results of Stage 2 and Stage 3 were generally comparable. However, estimates at Stage 3 were consistently and significantly lower than those of Stage 0, indicating the influence of ontogenetic adjustment. Interestingly, variations of NUM across different development stages exhibited an entirely opposite trend compared with the other variables.

Our modeling results clearly demonstrate that climatic variables are the primary drivers of stand structural dynamics, volume yield, and carbon accumulation in natural Larix gmelinii forests. As expected, for a given MAI, increasing SCI significantly enhanced Dg, BAS, VOL, and CAR across all developmental stages. This finding aligns with previous studies, which have shown that favorable climatic and edaphic conditions promote higher forest productivity and biomass accumulation in boreal ecosystems (Wang et al., 2011; Dong et al., 2024). Furthermore, the necessity of integrating climatic factors into empirical growth and yield models has been repeatedly emphasized (Lei et al., 2016; Trasobares et al., 2016; Zou et al., 2025), given that temperature, precipitation, and the corresponding drought indices strongly influence productivity trends and their spatiotemporal variability.

Interestingly, our simulations revealed that, for a fixed SCI, higher assumed MAI values (representing more humid conditions) were associated with lower Dg, BAS, VOL, and CAR. This counterintuitive result suggests that under stable stand structural conditions, a more humid climate (i.e., higher MAI) does not necessarily lead to greater biomass or carbon storage. Possible explanations include reduced solar radiation, lower temperature-driven growth rates and nutrient leaching from soils under excessively moist conditions (Mirabel et al, 2023), which collectively limit forest productivity. It reinforces that management goals focused solely on increasing increment rates without considering climatic or structural constraints may compromise long-term sequestration efficiency.

A notable strength of this study is its explicit accounting for the interaction between developmental stage and climate. The results reveal stage-dependent divergences in the trajectories of stand growth and carbon accumulation. Specifically, the predicted Dg, BAS, VOL, and CAR at Stage 1 were significantly higher than those of the stage-independent baseline (Stage 0), whereas Stage 3 exhibited a pronounced decline. These patterns reflect ontogenetic shifts in resource allocation, physiological efficiency, and self-thinning intensity throughout stand development (McMillan et al, 2008). Early-stage stands typically display higher responsiveness to favorable climatic conditions, while mature stands tend to experience physiological saturation and competitive constraints (Del Río et al, 2014).

The observed interactive effects suggest that the developmental stage acts as a critical mediator of climatic sensitivity. Similar findings have been reported in boreal and temperate forests, where stand age and structure strongly modify climate–growth relationships (McMillan et al, 2008; Wang et al, 2011). For example, Zhu et al. (2019) found that the climate sensitivity of net primary productivity (NPP) was greatest during dense stem phases following regeneration, illustrating that structural dynamics influence climate responsiveness. Our results thus confirm that growth and carbon sequestration in natural L. gmelinii are not static processes but evolve dynamically through developmental transitions.

The climate-stage interaction captured in our model offers valuable insights for optimizing forest management in line with China’s dual-carbon goals. The model indicates that younger and mid-developmental stages (Stages 1–2) exhibit higher growth and carbon accumulation potential compared with older stages (Stage 3), suggesting that stage-sensitive management interventions (e.g., density regulation) could maximize both timber yield and carbon sequestration. Conversely, older stands (Stage 3) demonstrated reduced growth and storage potential, implying the necessity of timely regeneration or structural rejuvenation to maintain sink strength. These findings align with evidence that forest carbon sequestration efficiency declines with stand age and that active management can extend the duration of high sequestration phases (Keenan, 2015; Li et al, 2025). In addition, partitioning total stand volume into assortments (large, medium, small, short, fuelwood, bark) and biomass carbon into different components (leaf, branch, stem and root) further supports rational planning of subsequent utilization and potential economic returns (Pukkala, 2018; Lin et al, 2024). This framework provides a foundation for adaptive, stage-sensitive decision-making that balances ecological and economic objectives in the context of climate change.

Despite the progress made, several limitations warrant attention. First, the model currently uses prescribed MAI values, which may not fully capture the effects of spatial and temporal variability in moisture conditions on stand growth and carbon accumulation. Second, disturbance factors (e.g., pests, fires) remain implicit, although such processes are expected to intensify under future climate change (Hu et al., 2012; Seidl et al., 2017). Third, while the framework is calibrated for L. gmelinii in Northeast China, recalibration for species- and region-specific applications is essential before its broad application. Future work should focus on integrating process-based submodules to capture physiological mechanisms and disturbance regimes, as well as employing remote-sensing products for large-scale validation. The development of coupled empirical–mechanistic frameworks will improve predictive accuracy and support more robust scenario analyses. Additionally, incorporating adaptive silvicultural options, such as variable retention thinning and assisted migration, could further enhance the model’s applicability and management decision-making.