Climate Change Effects On Height-Diameter Allometric Relationship Vary With Tree Species And Size For Larch Plantations In Northern And Northeastern China

: 10 Background: Tree height-diameter relationship is very important in forest investigation, understanding 11 forest ecosystem structure and estimating carbon storage. Climate change may modify the relationship. 12 However, our understanding of the effects of climate change on height-diameter allometric growth is still 13 limited at large scale. Methods: In this study, we explore how the climate change effects on height-diameter allometric 15 relationship vary with tree species and size for larch plantations in northern and northeastern China. 16 Based on the repeated measurement data of 535 plots from the 6th to 8th national forest inventory of 17 China, climate-sensitive tree height-diameter models of Larix plantations in north and northeast China 18 were developed by two-level nonlinear mixed effect (NLME) method. The final model was used to 19 analyze the height-diameter relationship of different Larch species under RCP2.6, RCP 4.5, and RCP8.5 20 climate change scenarios from 2010 to 2100. 21 Results: The values of 𝑅 𝑎𝑑𝑗2 (adjusted coefficient of determination), MAE(mean absolute error) and 22 RMSE(root mean squared error) of the NLME models for calibration data were 0.92, 0.76m and 1.06m, 23 respectively. The inclusion of climate variables MAT ( Mean annual temperature ), CMD ( Hargreaves 24 climatic moisture deficit ) with random effects was able to increase 𝑅 𝑎𝑑𝑗2 by 19.5% and reduce the AIC 25 (Akaike’s information criterion), MAE and RMSE by 22.2%, 44.5% Conclusion: According to the climate sensitivity, tree species could be classified as group I( L. gmelinii , 30 L. pincipis-rupprechtii and the unidentified species group) with large 𝛥𝐻 (from -4.77% to 18.17%) and 31 group II ( L. kaempferi and L. olgensis ) with small 𝛥𝐻 (from -6.37% to 9.4%).Large trees were more 32 sensitive to climate change than small trees. 33 Two-level climate-sensitive NLME model was developed for larch planatations in north and northeast China in this showed biological and statistical reasonability. MAT, CMD, the dominant climatic factor in modulating height-diameter allometry of larch plantations. Model simulation showed

principis-rupprechtii. In addition, there were trees are not identified to specific species which were 97 recorded as larch. According to the protocol of NFI, heights of 3-5 medium trees were measured in each 98 plot. In total, 7304 pairs of H-D measurements in 535 plots were obtained across seven Provinces. Data 99 were split into two parts for model calibration and validation by the following method: each plot was 100 randomly allocated to a number between 1 and 535, and plots whose number were less than 20th 101 percentile of all plots were assigned as validation data (1609 pairs in 107 plots) and the rest were fitting 102 data (5695 pairs of H-D measurements in 428 plots). Summary statistics of tree and stand variables can 103 be found in Table 1. The scatter plot can be found in Figure 1.   Climate data 124 The current climatic data for model calibration were downloaded from ClimateAP, which is an 125 application for dynamic local downscaling of historical and future climate data in Asia Pacific (Wang, 126 Wang et al. 2017). Seasonal and annual climate variables (averaged from 1980 to 2010) for a plot were 127 produced based on latitude, longitude, and elevation (Table 2). 128  Where is the total tree height (m), is the diameter at breast height (cm), BAL is the sum of basal 156 area larger than a subject tree, 0 , 1 , 0 , 1 and are model parameters, which have their own 157 biological characteristics, and is random error. 158 To evaluate the differences in height-diameter allometry among larch species, dummy variables S m 159 were created: (1) 1 = 1 denotes the L. gmelinii. and 0 the rest of cases; (2) 2 = 1 denotes the L. 160 olgensis and 0 the rest of cases; (3) 3 = 1 denotes the L. principis-rupprechtii. and 0 the rest of cases; 161 (4) 4 = 1 denotes the L. kaempferi; and (5) the category which can not be identified was represent by 162 1 = 2 = 3 = 4 = 0 as the reference. 163 Therefore, the model could be written as: 164 Where is the fixed-effect parameter vector, was dummy variable denoting tree species, and 166 other variables are defined as above. 167

Nonlinear mixed-effects climate-sensitive H-D model 168
To quantify the climatic effects on the H-D allometry, the selected climate variables were added into 169 the model by reparameterization for parameters in basic H-D model, and it could be written as: 170 Where was the climate variable vector selected by PCA and correlation analysis, and other 172 variables were the same as mentioned above.
Where and is the ℎ individual tree height nested within ℎ plot in the ℎ province, 180 and is the province-and plot-level random effects, and is the random error. Other variables 181 were the same as mentioned above. 182 The estimated random effect parameter were calculated as follows: 183 Where ̂ is the estimated prediction vector for random parameters, ̂ is the estimated q × q 185 variance-covariance matrix for among-unit variability, where q is the number of random effects 186 parameters in the model, ̂ is the estimated k × k variance-covariance matrix for within-unit 187 variability, ̂ is the partial derivatives matrix with respect to the random parameters and is the 188 residual vector determined by the difference between the observed and predicted heights using model 189 which only has fixed effects. 190 To account for the within-unit heteroscedasticity and autocorrelation in variance-covariance matrix 191 ( ), the variance-covariance matrix was determined as: 192 Where 2 is the value of residual variance of the estimated model, is a diagonal matrix 194 explaining the variance of within unit heteroscedasticity, Г is a diagonal matrix accounting for within 195 tree autocorrelation structure of errors, and AR(1) was used to reflect the within-tree autocorrelation 196 structure of errors for matrix Г . To reduce the heterogeneity in variance, the variance power equation 197 was determined as: 198 where, ̂ is the estimated height of ℎ tree nested in ℎ plot in ℎ province using fixed part 200 of the mixed-effects model; is the parameter to be estimated; and 2 is the same as defined in Eq.6. 201 Parameters in NLME models were estimated by restricted maximum likelihood implemented with 'nlme' 202 package in R software (Pinheiro, Bates et al. 2013). 203 When a new subject is available(for example, the one in the validation set), the model needs to be 204 Where n is the number of observations,̂ is ℎ tree estimated height of nested in ℎ plot nested 216 in ℎ province ; is the ℎ tree observed height nested in ℎ plot nested in ℎ province. ̅ 217 is the observed mean height for all data, , , are the total number of the province, the plots 218 nested in ℎ province, ℎ trees nested in the ℎ plot nested in ℎ province, is the number of 219 model parameters; and is the log-likelihood. 220 predicted under future and current climate scenarios, respectively. 233

Results 234
Selected climate variables 235 Three principal components described 95.27% of the variability of the climate data (Table 3) CMD were selected for further reparameterization using NLME. Summary statistics of MAT and CMD 244 can be found in Table 5.

Model comparison and evaluation 270
The fitting and validation results of the models are shown in Table 6. The base model (Eq.13) described 271 76% part of the variations in the height-diameter relationship when fitted the training data( 2 =0.76). 272 When tree species dummy variable and province-specific, plot-specific random effects were included in 273 the base model, the climate variables contributed significantly to the variance of tree heights and 2 274 increased from 0.77 to 0.92 (Table 6). Training dataset showed similar results, and the inclusion of 275 climate variables (Eq. 14) resulted in the increase 2 by 3% and the reduce of AIC by 3%. Mixed 276 effect model (Eq. 15) also removed the heteroscedasticity of residuals (Figure 2).  Results showed different effects of climate variables on parameters a 0 and b 0 , denoting the maximum 288 and relative change of tree height with diameter (Table 6). Parameter a 1 was significantly negative 289 indicating the increasing BAL will reduce the maximum height. The coefficient a 2 of MAT for parameter 290 a was significantly positive which means that the rising MAT will increase the maximum tree height. 291 This was also shown in Figure 3 where all H-D curves of different species became steeper under RCP2.6 292 and RCP 4.5 from 2010 to 2070. However, parameter b 2 was negative indicating that the rising MAT will 293 lower the tree height with the same diameter and there is a threshold for the effect of temperature on H-294 D relationship of larch species. Parameter a 3 was significantly negative indicating the decreasing 295 precipitation will reduce the maximum height. Both CMD and BAL showed marginal effects on H-D 296 relationship since b 1 and b 3 were nearly zero. 297 Table 6 showed that all parameters of tree species dummy variables f 1~f4 were positive, but f 1 andf 4 298 were not significant, indicating that L.olgensis and L.Kaempfer had significant difference with 299 unidentified group, which was also illustrated in Figure. 3. Coefficient f 2 was the largest indicating that 300 the maximum of tree height is the largest for L.olgensis. 301 MAT increases with the time and the temperature under RCP8.5 is largest followed by RCP4.5 and 302 RCP2.6. The precipitation under RCP8.5 is smallest and has the steepest slope followed by RCP4.5 and 303 RCP2.6. Figure 4 showed the H-D curve of larch species under climate scenarios RCP2.6, RCP4.5 and 304 RCP8.5. Generally, tree species can be obviously classified as two groups in terms of H, which are 305 group I (L. gmelinii group, L. pincipis-rupprechtii and the unidentified larch species) and group II (L. kaempferi under RCP2.6 and RCP8.5, and the sensivity was larger under RCP8.5 than that under RCP2.6. 320 However, the sensitivity was ranked as L. gmelinii > the unidentified species group > L. pincipis-321 rupprechtii > L. kaempferi > L. olgensis under RCP4.5. 322 2020). In this study, using mixed-effects model and including climate variables was able to increase 2 338 by 19.5% and reduce the AIC, MAE and RMSE by 22.2%, 44.5% and 41.8% for fitting set, respectively. 339 The residual heterogeneity was also reduced. Owing to the correlation among tree height-diameter II (L. kaempferi and L. olgensis). They showed strong ( H from -4.77% to 18.17%) and weak ( H from 381 -6.37% to 9.40%) response to future climate change, respectively. Under warmer and drier climatic 382 conditions, L.kaempferi and L.ogensis will grow thicker and shorter than the rest of the tree species group, 383 and their Hs are lower than those of group I for a given tree diameter. This may due to these two tree 384 species are moisture loving species (Wang et al., 1992). Under drought stress, the hydraulic conductivity 385 of the xylem of the trunk suffers irreversible loss. Therefore, the lack of water during the growing season 386 allows to allocate more resources for the growth of diameter (Ryan and Yoder 1997). Compared with 387 group II, group I is more resistant. As the temperature increases, more resources will be allocated to the 388 growth of the height than the diameter, thus trees would be higher. Previous studies also supported this Ethics approval and consent to participate 416 The authors declare that the study was not conducted on endangered, vulnerable or threatened species. 417 The authors declare that they obtained the informed consent from human participants involved in this 418 study. 419 Availability of data and material 424 The datasets generated during and/or analysed during the current study are available from the 425 corresponding author on reasonable request. 426