Mean Annual Wood Density Variations of Larix gmelinii (Rupr.), Quercus mongolica Fisch. ex Ledeb., and Pinus tabulaeformis Carr. at Two Di ﬀ erent Stem Heights

: Forests are a large carbon sink with an additional substitution e ﬀ ect in the merchantable timber compartment of harvested trees, where carbon stored within the same volume of wood varies depending on wood density. Here, we investigated mean annual air-dry wood density variations depending on cambial age, annual radial increment, and two di ﬀ erent stem heights of Larix gmelinii (Rupr.), Quercus mongolica Fisch. ex. Ledeb., and Pinus tabulaeformis Carr. from a ﬁrst climatic region (Mulan Forest) and exclusively of P. tabulaeformis from a second climatic region (Zhongtiaoshan Forest) in the temperate zone of China. We applied linear mixed-e ﬀ ects models with partly transformed variables and estimated marginal means for pairwise comparisons. Results showed that mean wood density was not signiﬁcantly di ﬀ erent between L. gmelinii (0.626 g cm − 3 ) and Q. mongolica (0.596 g cm − 3 ), but signiﬁcantly di ﬀ erent between P. tabulaeformis from the two di ﬀ erent climatic regions (0.445 g cm − 3 in Mulan Forest and 0.521 g cm − 3 in Zhongtiaoshan Forest). Mean annual wood density within trees except for P. tabulaeformis from Mulan Forest was initially increasing until an intermediate cambial age, after which it decreased again to lower values. These ﬁndings showed that tree age had to be considered in assessing carbon sequestration in wood. It also could play an important role in decision making for forest management in Mulan Forest and show the beneﬁt of the wood properties and carbon storage potential of the faster growing L. gmelinii compared to Q. mongolica . Furthermore, these ﬁndings gave an indication that intermediate old forest stands for some tree species accumulated more carbon per year within their woody biomass than young stands or old growth forests. Our results may have an impact on the planning of rotation lengths and of tree species composition for forest stands in Mulan Forest and Zhongtiaoshan Forest.


Introduction
Forest ecosystems are some of the largest sinks of terrestrial carbon on the globe [1]; they absorb billions of tons of atmospheric CO 2 yearly and contribute to climate change mitigation [2]. Within a forest, carbon is stored mainly in the living aboveground biomass and in the organic matter of the mineral soil [3]. The living aboveground biomass can be divided into several compartments, where the merchantable timber compartment has an additional carbon storage function when accounting for Boden et al. [36] studied the resolution abilities of HF densitometry and showed that, depending on the geometrical dimensions of the micro-electrode measuring system, it is possible to distinguish variations in wood density up to a resolution of 1 µm.
Here, we investigated mean annual air-dry wood density (d) variations depending on cambial age (CA), annual radial increment (ARI), and stem heights (SH) at 1.3 m (breast height (BH)) and between 10 m and 11 m (upper height (UH)) of L. gmelinii and Q. mongolica and of P. tabulaeformis from two different climatic regions. Air-dry wood density was determined by dividing the mass of air-dry wood by the volume of air-dry wood in the same state (other than for example basic wood density, where the mass of air-dry wood is divided by the volume of water saturated wood). We hypothesized that: 1.
Averaged over both stem heights and all observed CA, d was highest for Q. mongolica and lowest for P. tabulaeformis.

2.
For P. tabulaeformis, d had no significant difference between the two different climatic regions.

3.
With increasing CA, d increased for the two coniferous tree species L. gmelinii and P. tabulaeformis.

4.
For Q. mongolica, we expected d to decrease with increasing CA.

5.
For the effect of SH, we expected to observe significantly higher d at BH compared to at UH at the same CA for Q. mongolica and significantly lower d for the coniferous tree species.

Research Area
The research area of L. gmelinii, Q. mongolica, and P. tabulaeformis was located in Mulan Forest of Hebei Province (referred to as Mulan) and for the second location of P. tabulaeformis in Zhongtiaoshan Forest of Shanxi Province (referred to as Zhongtiaoshan) in Northeast China ( Figure 1). Both locations are dominated by temperate coniferous forests, which consist mainly of P. tabulaeformis at an elevation between 950 m and 1650 m above sea level. The climate at both locations can be classified as a temperate climate, which is influenced by the East Asian Summer Monsoon (EASM) with distinct hot and semi-humid summers and cold and dry winters [39]. Annual precipitation in Mulan ranges from 240 mm to 680 mm. Due to the EASM, about 70% of the total precipitation occurs in the wet season between May and July. A mean annual air temperature of 5.1 • C, ranging from −12.8 • C to 21.0 • C, characterizes this region (meteorological station of Weichang County, 844 m a.s.l., Hebei Province, 1951-2013). Annual precipitation in Zhongtiaoshan ranges from 331 mm to 885 mm, also with a strong influence of EASM. Mean annual temperature is 11.5 • C, ranging from −4.8 • C to 24.4 • C (aggregated climate data from Qinshui meteorological station, 887 m a.s.l., Yangcheng meteorological station, 661 m a.s.l., Gaoping meteorological station, 837 m a.s.l., and Lingchuan meteorological station, 1313 m a.s.l.). According to the Köppen climate classification system, the climate in Mulan is classified as Monsoon-influenced warm-summer humid continental climate (Dwb) and the climate in Zhongtiaoshan as Monsoon-influenced hot-summer humid continental climate (Dwa) [40].

Collection of Cross-Sections
Trees selected for stem disc collection were taken exclusively from the predominant and dominant canopy layer [41], as these trees experienced less influence from competitors on their growth, and their survival rate was higher compared to dominated individuals in a stand. In addition, sample trees were selected from closed single-species stands, which were not affected by gaps, skid trails, forest roads, and admixed tree species. Furthermore, selected sample trees had healthy and evenly shaped crowns and no visible injuries, damages, or infections (see Table 1 for additional information of sample trees). Moreover, selected forest stands were as old as possible to maximize the observation range of cambial age. Cross-sections were taken from felled trees at 1.3 m (BH) and between 10 m and 11 m stem height (UH), respectively. The collected cross-sections were then cut along one radius to facilitate a single controlled crack towards the pith while air drying for several weeks. Sampling took place during the growing season of 2014 in Mulan and during the growing season of 2018 in Zhongtiaoshan.

Density Measurement
Of each cross-section, four radii towards the four cardinal points were cut out into bars and glued on a glass plate. However, if a wood deformation due to an ingrown branch or injury was located directly on a radius towards a cardinal point and caused therefore a variation in tree ring structure and wood density that was not due to annual growth patterns, then the closest alternative radius without any deformation was measured instead.
For surface preparation of the wood bars, a single-point diamond fly cutter equipped with air bearings was used [42]. Here, a diamond blade that moved at a constant rotation speed along a circular track produced smooth surfaces by cutting small chips of the cell wall off the wood bars. Wassenberg et al. [43] showed that the use of the single-point diamond fly cutter with air bearings for wood surface preparation was superior compared to the preparation with a core microtome and sanding with P600 grit.
In this study, the HF densitometry was used to measure relative intra-annual wood density variations [32]. This method utilizes the dielectric properties of wood, which are closely related to wood density. Measurements took place via a micro-electrode system, where a slit-shaped probe had surface contact with a wood bar. The probe consisted of a high-frequency transmitting electrode and a receiving electrode, which were both hermetically separated by a thin metallic foil and had an integration width of 175 µm. The induced electromagnetic stray field propagated through the wood bar, and the detected voltage was higher in denser wood due to its higher local dielectric constant and was linked to the different ratios of cell wall to air-filled cell lumen. If temperature, moisture content, orientation of wood-fibers, and excitation frequency of the electric field are known, then the dielectric constant is linearly related to wood density [37,44].
The wood bars were mounted on a motorized linear travel stage that moved along a cross-sectional radius. The slit of the probe was oriented perpendicular to the measured radius. For HF densitometry, it is crucial that the surface preparation of the wood bars fulfills a high standard, ensuring that no air gaps at any point during measurement are present and that the probe has always full contact with the wood bar. Otherwise, the gap of air between the probe and the wood bar will lead to a decrease in the output voltage and therefore implies lower density [32,43].
Both the single-point diamond fly cutter equipped with air bearings and a prototype of the HF densitometer were developed and set up at the Chair of Forest Growth of the Albert-Ludwigs-University Freiburg, Germany. The utilized software was developed at the Chair of Forest Growth of the Albert-Ludwigs-University Freiburg, Germany, which enabled superimposition of wood density profiles to tree ring series.

Data Processing
HF measurement of all cross-sectional radii took place over several days. To normalize all measurements to a date where a calibration function was calculated, a reference radius was measured at each day of measurement. The mean density of the total daily reference radius (x R i ) was then calculated and divided by the mean density of the total reference radius from the day when the calibration function was calculated (x R cal ). Each measured relative density value was multiplied by the corresponding daily quotient: where d norm i is the normalized relative volumetric density value at the day of HF measurement i, x R i is the mean relative volumetric density value of the total reference radius at the day of HF measurement i, x R cal is the mean relative volumetric density value of the total reference radius at the day of HF measurement when a calibration function was calculated, and d i is the relative volumetric density value at the day of HF measurement i. The calculated daily quotients ranged from 0.950 to 1.027.
Calibration functions were repeatedly determined to convert measured relative volumetric density from HF densitometry to true volumetric mass density. These calibration functions were based on the linear relationship between the measured relative volumetric density of 6 tropical wood specimens and their true mass density, which was known by determining the weight and volume of the air-dried wood specimen. By use of the calibration functions, the values of voltage can be calculated into values of density: where y is the density of air-dry wood in g cm −3 , a is the parameter for the slope of the linear relationship, b is the parameter for the intercept at HF output of 0, and x is the measured value of relative density from HF densitometry in volts.
For each density profile, the inter-year transient was removed, as it was an artifact from the integration width of the slit-shaped probe. Specifically, within the data points of a single year, all data points up to the first local minimum and after the last local maximum were removed (see Figure 2 for an example). Of the cross-sectional radii, mean annual density was calculated as the arithmetic mean of all measured values. d was only analyzed from the ring of CA 2, as the innermost tree ring of the cut wood bars was in many cases only fragmented due to the width of the saw, and therefore, no true mean values of annual gravimetric density could be obtained. Furthermore, sampling of stem discs took place during the growing season, and the last tree ring of the current calendar year was also not analyzed, as it could never be assured that wood formation processes of the current year had already finished.

Data Analysis
All statistical methods and data analysis were computed in the R programming environment using the Graphical User Interface R Studio [45,46]. d was analyzed by using linear mixed-effects models (LMEM), which are especially useful for analysis of non-independent data from repeated measurements with sporadically missing values. The LMEMs were calculated for species subsets using the lmer function of the package lme4 [47]. As scatter plots of mean annual air-dry wood density for each tree species plotted over CA and ARI revealed partly parabolic patterns and partly uneven variability (Figure 3), LMEMs were applied for subsets of species and location, and variables were partly transformed. The models were generally formulated as: where D i,j,k,l is the response variable mean annual air-dry wood density (d) of the corresponding species and location (L. gmelinii, Q. mongolica, and P. tabulaeformis from Mulan and P. tabulaeformis from Zhongtiaoshan), µ is the general experimental mean, CA i is the continuous fixed-effect CA of cambial age i, ARI j is the continuous fixed-effect ARI of annual radial increment j, SH k is the categorical fixed-effect SH of stem height k (BH and UH), and (CA × ARI × SH) i,j,k are the interaction terms, I l is the random-effect I of the individual tree l, and ε i,j,k,l is the residual error term. To obtain more normal residuals, for L. gmelinii, D i,j,k,l was applied with a power transformation using the exponent of 0.5, and CA i was formulated as a second order polynomial term of CA; for Q. mongolica, CA i and ARI j were formulated as second order polynomial functions of CA and ARI; for P. tabulaeformis from Mulan, D i,j,k,l was applied with a power transformation using the exponent of 0.5; and for P. tabulaeformis from Zhongtiaoshan, D i,j,k,l was applied with a power transformation using the exponent of 0.75, and CA i was formulated as a second order polynomial term of CA. Diagnostic plots of the residuals of each LMEM for species subsets are shown in Appendix A Figures A1-A4 and were satisfying. To calculate ANOVA (analysis of variance) tables with p-values and F-statistics, the lmerTest package was used, where the anova and summary functions for objects of the lmerMod class were overloaded to provide p-values for fixed-effects [48].
Post-hoc tests between subsets were performed as pairwise comparisons using the function emmeans of the package emmeans [49]. To avoid higher ranking interactions and therefore potentially misleading results, the model was formulated for the whole and untransformed dataset with the subset (S) of species and location as fixed-effect and cambial age (CA), stem height (SH) and individual tree (I) as random-effects. Diagnostic plots of the residuals are shown in Appendix B and were satisfying. For visualization of interaction effects with continuous predictor variables of the LMEMs, plots were produced based on effect displays obtained by the function Effects in the package effects [50].

Effect of Tree Species
A significant effect of subsets (p < 0.001) was caused by significantly lower d for P. tabulaeformis from Mulan (0.445 g cm −3 ) compared to all other tree species (L. gmelinii, 0.626 g cm −3 , p < 0.001; Q. mongolica, 0.596 g cm −3 , p < 0.001; P. tabulaeformis from Zhongtiaoshan, 0.521 g cm −3 , p = 0.050). d of P. tabulaeformis from Zhongtiaoshan was also significantly lower than d of L. gmelinii (p = 0.005), but no significant difference of d between Q. mongolica and L. gmelinii and of d between Q. mongolica and P. tabulaeformis from Zhongtiaoshan was observed ( Figure 4).

Effect of Stem Height
According to the LMEMs for each subset, a significant effect for SH was only found for L. gmelinii and P. tabulaeformis from Mulan (Table 2). For L. gmelinii, d at BH was on average 0.656 g cm −3 and at UH on average 0.625 g cm −3 . For P. tabulaeformis from Mulan, d at BH was on average 0.463 g cm −3 and at UH on average 0.436 g cm −3 .

Effect of Cambial Age
The fixed-effect of CA was significant for all subsets ( Table 2). As can be seen from the effect plot of d over CA for all subsets, no clear linear increasing or decreasing age trends were visible, but partially parabolic patterns were present in the data ( Figure 5). The interaction of CA with SH for subsets was significant for all groups, but Q. mongolica ( Table 2). The effect plots for the predictor variables CA and SH of subset models showed similar patterns for both SH of L. gmelinii and P. tabulaeformis from Mulan (Figure 6a,b). The significant interaction effect of CA and SH of P. tabulaeformis from Zhongtiaoshan showed similar density values for the first 40 to 50 years. However, the fitted model indicated a convex and overall increasing pattern for d at UH, while at BH d, was overall decreasing in a concave pattern (Figure 6c).

Effect of Annual Radial Increment
ARI had a significant effect on d for all subsets, but P. tabulaeformis from Zhongtiaoshan ( Table 2). The effect plot of d over ARI revealed an increasing linear pattern for d with increasing ARI for L. gmelinii and a decreasing linear pattern for P. tabulaeformis from Mulan. For Q. mongolica, increasing d with increasing ARI was observed for ARI below 3 mm, after which, it leveled off and decreased again (Figure 7).

Effect of Cambial Age, Annual Radial Increment, and Stem Height
The three-way interaction effect of CA, ARI, and SH was significant for all subsets ( Table 2). Effect plots of the predictor variables CA and SH for multiple levels of ARI showed for L. gmelinii that the parabolic pattern at BH was diminishing with increasing ARI and that a close to linear increase was present. At UH, the curve of d over CA at an ARI of 1 mm had an overall decreasing shape between 15 years and 72 years. At ARI of 2 mm and 3 mm, a parabolic shape was present and shifting towards higher d with increasing ARI. At ARI of 4 mm and 5 mm, an increase in d for observed CA was present (Figure 8a-f). For Q. mongolica, the corresponding effect plots revealed a different development. At BH, the parabolic pattern of d with an increase in the first 20 to 25 years remained over all levels of ARI, but was increasing on average. At UH, the curve at ARI of 1 mm was convex, whereas for ARI of 2 mm, the shape for the observed range was concave with a steep slope (Figure 9a-c). Effect plots for P. tabulaeformis from Mulan showed a linear decrease of d with increasing CA at BH for tree rings with ARI of 1 mm and almost constant values of d at UH. With an ARI of 3 mm, the relationship changed to linearly increasing for both SH, which became stronger for higher values of ARI and which was also always more pronounced at BH than at UH (Figure 10a-e). For P. tabulaeformis from Zhongtiaoshan, the corresponding effect plots revealed a concave parabola at BH that was especially visible at ARI of 2 mm and 3 mm. At UH, modelled d shifted from a convex pattern for ARI below 3 mm towards a more linear increase at ARI of 3 mm and 4 mm (Figure 11a-e).

Discussion
In this study, we used linear mixed-effects models and pairwise comparisons of estimated marginal means to analyze variations of mean annual air-dry wood density (d) of L. gmelinii, Q. mongolica, and P. tabulaeformis from Mulan and P. tabulaeformis from Zhongtiaoshan over CA and ARI and between breast height (BH) (1.3 m) and upper stem height (UH) (between 10 m and 11 m).
One of our main findings was that d was significantly higher for L. gmelinii compared to P. tabulaeformis from both locations and that Q. mongolica was only significantly higher than P. tabulaeformis from Mulan, but not significantly different from L. gmelinii and P. tabulaeformis from Zhongtiaoshan ( Figure 4), which was not fully in line with our first hypothesis. We formulated our first hypothesis according to Zanne et al. [6], where average wood density for Q. mongolica was between 0.603 g cm −3 and 0.660 g cm −3 , average wood density for L. gmelinii between 0.528 g cm −3 and 0.599 g cm −3 , and average wood density of P. tabulaeformis 0.360 g cm −3 . However, the values of d for L. gmelinii found in this study were higher than most values from the literature [6,51], while those values of d for Q. mongolica were lower than most values from the literature [6,24]. Wassenberg et al. [52] showed for multiple tree species that, as wood density values of a single tree can vary within the stem, estimated mean values also vary based on the sampling strategy. It is therefore possible that under a different sampling strategy, mean wood density values would show higher similarity to cited values.
Other than hypothesized, the means of d were significantly higher for P. tabulaeformis from Zhongtiaoshan than those of P. tabulaeformis from Mulan. Figure 3g,h shows that different ARI for the different locations could not be the cause as d remained systematically higher at Zhongtiaoshan over the whole observed range of ARI. Furthermore, the longer observation period of CA in Mulan of 100 years was not an explanation for the observed differences of d as also for this predictor, d remained relatively stable ( Figure 5). A possible explanation could be differences in climatic conditions, as meteorological data indicated a warmer and moister climate at Zhongtiaoshan, which in turn may had led to adaptions to different site conditions. In fact, when looking at average tree height at 50 years of tree age, trees in Mulan were 15.9 m high and trees in Zhongtiaoshan 12.2 m high, indicating different site productivities. Gindl et al. [53] found that lignin content in secondary cell wall layers of terminal latewood tracheids of P. abies positively correlated with the air temperature of August and September. Rossi et al. [54] suggested that cell production of conifers in the early growing season allowed tracheid differentiation to be completed and cell wall lignification to be adequate. As the growing season in Zhongtiaoshan was longer and started earlier than in Mulan with higher monthly air temperature, a higher amount of lignification and therefore denser wood could be the result. Larson [55] found on Pinus elliottii Engelm. in the south-east United States that the percentage of latewood within a tree ring increased with higher precipitation in June and July. Additionally, the evolution of different provenances could also explain differences in d between the two locations, leading to small differences in wood anatomy and growth responses. Follow-up research with the consideration of climatic parameters and an analysis of growth response to climatic conditions could provide explanations of different wood density values for P. tabulaeformis. Furthermore, additional research on wood density depending on different site productivity classes could lead to explanations of our observed results.
Contradicting to our third hypothesis, we did not always observe obvious increasing patterns of d with increasing CA for the coniferous tree species ( Figure 5); only at UH of P. tabulaeformis from Zhongtiaoshan, an overall increase was present (Figure 6c). The expected reason was that with increasing CA, ARI decreased, while at the same time, the ratio of latewood to earlywood within a tree ring increased [56]. The majority of cells in the xylem of coniferous tree species are tracheids, which have a higher ratio of cell wall to lumen in latewood, causing a higher wood density compared to earlywood. However, for P. tabulaeformis from Mulan, decreasing d with increasing ARI could be observed, which was expected to be linked with increasing d over increasing CA (Figure 7). This finding did not correspond with the results of Xu [19], who observed an increase of wood density after the 10th tree ring for P. tabulaeformis. Possibly an increase in sample size and a follow up study could lead to clearer results.
An expected decrease of d for increasing CA of Q. mongolica could be observed only after an initial increase ( Figure 5). A Similar decrease of wood density with increasing CA was shown for other oak species such as Quercus suber L. [57] and Quercus petraea Liebl. [58]. As a ring-porous broadleaved tree species, Q. mongolica produces large diameter tracheas for water transportation in the beginning of the growing season. With narrower tree rings, the ratio of tracheas within the xylem increases, resulting in a decrease of wood density [56]. As expected, d increases with increasing ARI, which was expected to be linked with decreasing d over CA (Figure 7). The same relationship of increasing wood density with increasing tree ring width was found by Zhang et al. [59] for Q. petraea and Quercus robur L. This led to the conclusion that for the investigated trees, no clear linear relation between CA and ARI was present. Especially at UH with an average observation period of 28 years, long term growth may not be detectable. However, a parabolic pattern in specific gravity over CA was observed for the closely related tree species Quercus liaotungensis Koidz., which supported our results [60].
For the effect of SH, we expected to observe significantly higher d at BH compared to at UH at the same CA for Q. mongolica and significantly lower d for the coniferous tree species. For annual growth during the same calendar year, different studies for multiple tree species showed that no differences in ARI were present, explained solely by different stem heights [61,62], and that ARI was highly similar at different stem heights [63]. We assumed therefore that ARI at younger CA at UH was less than ARI at younger CA at BH. For tree rings with less ARI, we expected less d for Q. mongolica and higher d for both coniferous tree species, as explained already in relation to tree ring anatomy.
A significant effect of SH was only present for L. gmelinii and P. tabulaeformis from Mulan. Figure 3 indicates that d was lower at UH of both subsets. Repola [64] found for P. sylvestris in southern Finland also a decrease of wood density with increasing SH. Furthermore, the pattern of decreasing wood density with increasing SH was observed for Pinus massoniana Lamb. and Pinus koraiensis Siebold & Zucc. in China [16,65]. However, Jyske et al. [12] found for P. abies an increasing pattern of wood density with increasing SH in a given CA, which showed that the results in the literature were not consistent for the family Pinaceae.
A contrasting development of d with increasing CA was observed for the two different SH of P. tabulaeformis from Zhongtiaoshan; at BH, d overall decreased in a concave pattern, while at UH, d increased in a convex pattern (Figure 6c). Additionally, Figure 11 reveals that the initial convex pattern at UH changed towards a more linear increasing pattern for higher values of ARI. Molteberg and Høibø [66] found a similar convex pattern for basic wood density over cambial age in 28-year-old P. abies stems. Furthermore, Petty et al. [67] found for P. abies and Picea sitchensis (Bong.) Carr. an initial decrease of wood density followed by a gradual increase. However, for all curves at UH, data were present for relatively short observation periods, leading to some degree of uncertainty in the modelled values.

Conclusions
Our study showed that mean wood density was not significantly different between L. gmelinii and Q. mongolica. This finding could play an important role in decision making for forest management in Mulan and showed the benefit of the wood properties and carbon storage potential of the faster growing L. gmelinii compared to Q. mongolica, as indicated by overall narrower ARI for Q. mongolica and overall higher sample trees of L. gmelinii (Figure 3e,f, Table 1). Furthermore, we found that wood density within a tree was not monotonically increasing with cambial age, but for some subsets, levelling off eventually and even decreasing to lower annual values. This finding gave an indication that intermediate old forest stands possibly accumulate more carbon per year within their woody biomass than young stands or old growth forests, while at the same time, the time of maximum mean annual carbon sequestration would be delayed. This conclusion may have an impact on the planning of rotation lengths for forest stands in Mulan and Zhongtiaoshan. However, it also showed that tree age needs to be considered in assessing carbon sequestration in wood.
Furthermore, we stress that due to the small sample size of five trees per subset, the external validity of our results was not given, and the results were only valid for pre-dominant and dominant trees of the same tree species in the same geographical region, Mulan and Zhongtiaoshan. Nevertheless, our work served well as an exploratory study, indicating interesting and novel results and interactions of effects on wood density variations. We therefore see high benefit in follow-up research on wood density variations within stems of tree species in China with increased sample size and possibly additional levels of the effect of stem height. Furthermore, additional attention needs to be given to a climate analysis of growth reactions and the effect of site productivity and tree age on wood density.