Process-based modeling analyses of Sabina przewalskii growth response to climate factors around the northeastern Qaidam Basin

Sabina przewalskii is the longest living endemic tree species in the northeastern Tibetan Plateau, and has been widely employed in tree ring based climate research in China. However, most dendroclimatic reconstructions have been primarily based on empirical relationships between tree growth and climate factors identified by statistical assessment. To date, the physiological relationships between tree growth and their limiting climate factors have not been properly assessed. Here, we simulated the physiological response of Sabina przewalskii tree growth to major limiting climate factors based on the Vaganov-Shashkin (VS) model. The VS modeled results validated the relationships between tree ring and climate factors constructed by statistical models, both approaches suggesting that precipitation during the early growing season, especially in May and June, has significant effect on tree growth, while temperature mainly affects tree growth by warming-induced drought and by extending the growing season. Under current and projected climate scenarios, our modeling results predict an increase in radial growth of Sabina przewalskii around the Qaidam Basin, with the potential outcome that regional forests will increase their capacity to sequester carbon.

Tree ring indices (width, maximum latewood density and isotope etc.) play a crucial role in paleoclimatic research [1,2]. Several millennial climatic reconstructions at regional and hemisphere scale have been developed by this proxy [3,4]. Tree ring indices are useful in revealing historical climatic variation and in projecting future climate change [5]. One common approach in dendroclimatology is to approximate the relationship between tree growth and climate factors by a linear regression function [6,7]. However, as temperatures have increased in many parts of the Earth during the last century, a phenomenon called "divergence" between tree growth and climate factors has been reported from the Northern hemisphere circumpolar area and from mountain area in the mid-latitudes of the Northern Hemisphere. This phenomenon encompasses the divergence be-*Corresponding author (email: shaoxm@igsnrr.ac.cn) tween tree ring based reconstruction and measured data [8][9][10], the divergence of tree growth and their main limiting climate factors [11,12] and the divergent growth responses of trees within sites [13,14]. All of these divergences indicate that the statistically approximated linear relationship might not universally apply between tree growth and their limiting climate factors. Tree ring based climate reconstructions based on statistical models might over-or underestimate climate history, highlighting the need for a systematic physiological analysis of tree growthclimate relationships.
Physiological tree growth models are process-based models based on linear and non-linear relationships between tree growth and climate factors. Among such models, the Vaganov-Shashikin (VS) model is the most popular and widely used [15,16]. The VS model focuses mainly on the relationships between cambium activity and its limiting climate factors. The model is based on the following hypotheses: (1) External influences mainly affect the cambial zone, specifically, the linear growth rate of cambial cells.
(2) The external factors are temperature, light and soil moisture. (3) The principle of limiting factors is employed in the calculation of growth rates i.e. growth rate at a certain time of a season cannot be higher than that allowed by the most limiting factor. (4) The variations of growth rate in the cambial zone mainly pre-determine the anatomical characteristics of the tracheids being formed (e.g. radial diameter, cell wall thickness). (5) The model simulates variations in growth rate and tree ring structure resulting from current climatic changes (within the season). However, the model simulates the climatically induced tree ring width and structural variations only. Compared to other models, e.g. TREERING2000, that explicitly treats photosynthesis, respiration and the partitioning of assimilates [17], its applicability is limited. Despite such limitations, the VS model remains the most popular process based physiological model on account of its simplicity and its proven success in simulating coniferous tree growth in North America and Siberia [18,19]. In China, the VS model successfully simulates the growth of Pinus tabulaeformis from Helan Mts and Pinus armandii from Huashan [20][21][22], and was hence an obvious choice for the simulation of radial tree growth in our current study species.
Sabina przewalskii is an endemic conifer species on the northeastern Tibetan Plateau [23]. It is the most important forest component because of its long age, drought endurance, and wind resistance. Several millennia-long tree-ring width chronologies of Sabina przewalskii have been compiled [24][25][26][27] and employed in regional precipitation recon-struction. However, the tree growth climate relationships underlying these reconstructions were mainly assessed by statistical models, and no process based physiological validation has been undertaken to date. In this study, we simulated Sabina przewalskii growth and its main limiting climate factors by the physiological process model (VS model) to test (1) how trees physiologically respond to variations in their main limiting factors, (2) whether results are comparable to this relationship as derived by statistical models, (3) how this relationship has varied in the past and what it predicts for the future.

Tree ring data
Thousands of Sabina przewalskii tree cores from seven sites have been collected from Zongwulong Mt. and Shalike Mt. around the northeastern Qaidam basin ( Figure 1). Increment cores were taken from dominant and co-dominant trees which appeared healthy and which were isolated at or close to their upper limit of around 3500-4000 m a.s.l. [24]. All cores were processed following standard methods [28]. All series have been employed to reconstruct the millennial precipitation in the Delingha region. After comparing the relationship between the sites, we chose DLH3 (37°27′05′′N, 97°32′33′′E) as a typical site representative of all seven sites to simulate the regional physiological process of Sabina przewalskii to its main limiting climate factors.
DLH3 is the closest sample site to Delingha city, in which the meteorological station is located. 56 cores from 28 trees appeared healthy and isolated were sampled at a northwest facing site with slope around 20°, 50 cm soil layer at 3920 m a.s.l. Following precise crossdating, chronologies were developed using the ARSTAN program [24]. Nine pointer years were selected from 11 trees according to the annual precipitation and common wide rings and narrow rings. Wide rings developed in 1959, 1989, 1993 and 2002, while narrow rings were formed in 1957, 1966, 1979, 1995 and 1998. Cell numbers of each pointer year were counted by image analysis. First, the mean xylem cell numbers from five files were counted in each ring of every pointer year on each selected tree. The cell number of a given pointer year is the average of all standardized mean cell numbers of each tree in that year [28].

Climate data
The continental arid climate pronounced in our study area is primarily influenced by the Westerlies. The average monthly and daily temperatures and precipitation data from the meteorological station in Delingha (Figure 1), which holds climate records dating from 1955 (1956 in our analysis because some data from 1955 are missing) have been used in this study. The Mann-Kendall test was applied to test the significance of trends in annual temperature and precipitation. Because the elevation of the meteorological station (2981.5 m) is lower than that of the sample site (3920 m), a lapse rate (0.56°C/100 m) was applied to interpolate the measured temperature from Delingha station to our sample site. No interpolation was performed on the precipitation data.

VS model
In the model, the seasonal growth of a tree is assumed depended on three main factors: where GrE (t), GrT (t) and GrW (t) are the partial growth rates, calculated independently from solar irradiation, temperature and water content in soil.
(1) Calculation of solar radiation on tree growth. GrE (t) is calculated from eq. (2), where ϕ is latitude, θ is sun slope angle, and ϕ is day length.
Here the transmittivity of the atmosphere and the eccentricity of the Earth around the sun are ignored.
(2) Calculation of temperature on tree growth. In the model, tree growth starts when accumulated temperature is equal to T beg . A piecewise linear function is employed for calculation of temperature on tree growth. There is no growth when temperature is less than a given minimum temperature. The growth rate increases linearly up to the first temperature optimum (T opt1 ). Growth rate is stable between T opt1 and the second temperature optimum (T opt2 ) in-clusive. The growth rate decreases linearly between T opt2 and a limiting temperature (T max ), beyond which growth ceases.
(3) Calculation of soil moisture on tree growth. Similarly to GrT (t), GrW (t) is calculated by a piecewise function with four parameters W min , W opt1 , W opt2 and W max . Daily water in soil (d W ) is determined from a water balance equation: where f (P) is daily precipitation, Er is daily water loss in soil by transpiration and Q is runoff. f (P) is calculated as where P is the actual daily precipitation and k1 is the part of precipitation falling into soil. There are two ways to evaluate and define the parameters in the model. All parameters can be derived from experimental information prior to calculation. Alternatively, approximate starting values of parameters can be used in initial calculations and the simulated results are then compared with the observed tree ring chronologies, as the first step in an iterative model-improvement process. In this study, both approaches have been applied. Parameters such as the depth of root system, wilting moisture and maximum daily precipitation falling into soil were defined by observed information. The other parameters were determined using the iterative approach. Start values used here were the same as those in previously-published studies [21].

Climate condition in study area
The continental arid climate is pronounced in the study area, which is characterized by a typical inland aridity. The warmest and wettest month is July. Throughout the past 50 years, the highest recorded temperature in July was 20.5°C and the maximum precipitation was 36.1 mm. The main source of soil moisture is precipitation. According to climate data recorded at Delingha meteorological station (Figure 2), annual mean temperature and precipitation have increased significantly during the last half century (P<0.0001). Temperature increased mainly during autumn and winter (Jan., Feb., Sep., Nov. and Dec.), while precipitation increased mainly in July (P=0.0008), January (P=0.01) and June (P=0.09). These results suggest an upward warming and wetting trend in this region.

Statistical analysis
The standing and environment of Sabina przewalskii is similar at northeastern Qaidam Basin. The climate information recorded by trees sampled within this area is coincident also. According to correlations between chronologies after removing growth trend and retaining most of the low frequency signals (Table 1) [24], it is clearly seen that trees from these sites have been coherent at both low and high frequency (first difference) during the last thousand years. The highest correlation (0.827) is between DLH1 and DLH2, the lowest (0.598) is between DLH1 and WL4. Correlation relationships of these two pairs are significant at first difference with correlation coefficients 0.919, 0.696 respectively. DLH3 shows high correlation relationships with all other sites. Hence we used DLH3 to represent all the other sites in the following analysis.
Trees at the DLH3 site have high common variance. The mean sensitive is 0.4 and miss ring rate is 1.238%. The ring width chronology is 1599 a, AD404-2002. The statistical model shows that DLH3 ring width chronology correlates positively with precipitation of current May (R=0.38, P<0.01) and June (R=0.65, P<0.001) and negatively with temperature of current June (R= -0.39, P<0.01, Figure 3).

Process model analysis
(1) Model parameters. Tree ring model parameters used in our final simulation are listed in Table 2. Model parameter testing revealed that W opt1 , T min and T opt1 are the most sensitive  parameters to tree growth (Figure 4). The correlation between the observed and calculated values of tree ring width increases as W opt1 increases in the range 0.1<W opt1 < 0.18. This trend disappears for values of W opt1 >0.18. No obvious change is observed when other parameters are changed (e.g. T opt2 , W opt2 ).
(2) Model behavior: inter-annual and intra-annual simulations. We first assessed the ability of the VS model to simulate inter-annual tree ring formation. The observed and simulated annual tree ring indices of the analyzing period, 1956-2002, are plotted in Figure 5. In general, good agreement is observed at annual and at inter-annual scale. The correlation between observed and simulated series during the analysis period is 0.681 (P<0.001). Xylem cell numbers of pointer years also match well between actual and simulated data, with correlation 0.73 (P<0.01). Comparing model result and observed data of single trees, the highest correlation is 0.92(P<0.001) and the lowest is 0.50 (P< 0.01).
(3) Influence of climate factors on ring width formation. The difference between a wide ring year and narrow ring year is caused mainly by moisture limitation (Figure 6). A high annual radial growth (wide ring formed) results when the moisture is abundant and temperature is not high during the beginning of the growing season (May to June). On the contrary, small annual radial growth (narrow ring formed) results when moisture is reduced during the beginning of the growing season, although the growth rate due to temperature Figure 4 The sensitivity range of model parameters.  is high at the same time. According to the observed climate record, no large differences of annual precipitation (greater than 100 mm) exist between pointer years, excepting 1957 (during which the annual precipitation was 70 mm). The main difference between years results from the accumulation of precipitation from January to June, especially in April, May and June. In years when a wide ring is formed, the precipitation during April-June is around 150 mm, about 50% to 75% of the annual total precipitation. In years with narrow rings formed, the precipitation during April-June is less than 50 mm, a mere 28% of annual precipitation. Meanwhile, there is no difference in temperature during this period between all pointer years.
(4) Effects of growth start and end dates on cell numbers and tree ring widths. According to our model simulation (Figure 7), the start date of the growing season falls within the range of the 114th to the 218th day of each calendar year (corresponding to end of April to beginning of August). The end date of the growing season falls within the range of the 234th day to the 284th day (corresponding to end of August to beginning of October). The length of the growing season of Sabina przewalskii has expanded in recent years. In the last half of the analyzing period, the start date is advanced by 9 days, and the end date delayed by 6 days, compared to the first half part of the analyzing period. Correlation coefficient for the observed tree ring index series and the start date series is -0.611, which is significant at the 0.001 level. There is no significant correlation between the observed tree ring index series and the end date series.

Discussion
The VS model has been successfully employed to simulate, evaluate and interpret the tree growth-climate relationship of several conifer species for a variety of environmental conditions in the northern hemisphere [17,19,21,22]. Here, the VS model performed well in simulating the growth-climate relationship of Sabina przewalskii in a dry, high-altitude area. High correlations between observed and simulated values indicate that the VS model could simulate the radial growth of Sabina przewalskii well at both the annual radial growth and at the cell level. Hence, it can be used with confidence to evaluate the physiological relationship between Sabina przewalskii tree growth and its major limiting climate factors.
Physiologically, different species growing in different environments respond differently to climate factors, as does the same species growing in different environments. Consequently, parameters in the model are characterized by different species and sites. Our values of W opt1 , T min and T opt1 are lower than those used in the simulation of Pinus tabulaeformis from Helen Mt. [22]. (for which W opt1 =0.4, T min =5°C, T opt1 =20°C). Reasons for this difference might be related to the habitat of Sabina przewalskii. Sabina przewalskii grow mainly on sunny slopes of mountains surrounding Qaidam Basin, where more solar radiation is received than on other slopes and trees can start growth at lower temperature. Meanwhile, trees are expected to utilize all available moisture to grow in this arid area, especially given the high evaporation induced by enhanced solar radiation. On the other hand, parameters values used in the simulation of Sabina przewalskii are similar to those of some conifer species at boreal tree line (e.g. Yakutia, T min =4°C, 15°C<T opt1 <20°C) [17]. This finding could be explained by the high and cold climate condition of the Tibetan Plateau, which is similar to that found at high lati-tudes.
The climate conditions during the first part of the growing season greatly influence ring width formation at midand high latitudes [17,29,30]. Precipitation during the beginning of the growing season is crucial for radial growth of Sabina przewalskii. This is reflected in results from the statistical model, which reveal a positive correlation with precipitation during May and June, and a negative correlation with temperature in June. In the physiological model, wide ring formation and more xylem cell growth is associated with years of increased precipitation during May and June (e.g. 1989, 1993 and 2002). The effect of temperature is lower than that of precipitation on the radial growth of Sabina przewalskii in our study. Temperature mainly affects tree growth indirectly by warming-induced drought at the beginning of the growing season. Although temperature has increased during the last century in our study area and the growing season has extended, the influence of temperature on growth rate has not changed over this period.
The divergences phenomenon between tree growth and main climate factors is one of the most important research interests in dendroclimatology [8][9][10]. It has been reported not only at high latitudes in the northern hemisphere [9,10] but also in mountain areas of Eurasia [30,31]. The divergence between tree growth and growing season temperature is potentially driven by warming-induced moisture stress [10,13,14]. The crucial role of precipitation on the radial growth of Sabina przewalskii in Qaidam Basin has been demonstrated in both statistical models and a physiological model. However, according to our research, recent warming has not affected the tree growth climate relationship in recent decades.
During the past 50 years at Qaidam basin, significant temperature increases have occurred mainly in autumn and winter, while precipitation has increased primarily during summer and winter [32][33][34]. This change of hydrothermal condition is advantageous to tree growth. Comparing the growth rate in two parts of the analyzing period (Figure 8), the growth rate due to precipitation (GrW) in the second part is greater than that in the first part. No obvious change in growth rates due to temperature (GrT) has occurred between the two time intervals. The total growth rate (Gr) is accelerated in summer and is also affected by moisture condition, consistent with observations that the mean ring width in the last half of the analyzing period is wider than that of the first half. Recent studies have projected an ongoing warming and wetter trend in northeastern China [35]. The results presented here suggest that this trend has already begun in the Qaidam Basin. Under current and projected climate conditions there will be an increase in radial growth of Sabina przewalskii around the Qaidam Basin. An important potential outcome of this prediction is that regional forests may enhance growth and increase their capacity to sequester carbon.

Conclusion
Using the VS model, we simulated the physiological response of Sabina przewalskii tree growth to its major limiting climate factors. Results show that this model can successfully simulate the radial growth of Sabina przewalskii at both ring width and at cell level, and can reveal the effects of different climate factors on tree growth at different stages. Precipitation during the early growing season, especially during May and June, has significant effect on tree growth, while temperature mainly affects tree growth by warming-induced drought. No change has been found in tree growth climate relationships over the last half century. However, the growing season of Sabina przewalskii has been extended as a result of recent warming. The radial growth has also increased, because of higher precipitation at the beginning of the growing season during recent decades. Our modeling results also suggest future favorable growth conditions for Sabina przewalskii around the Qaidam Basin, leading to a potential increase of local carbon sequestration.