Climatic Gradients of Biomass and Net Primary Production of Mixed Picea-Abies Forests in Eurasia

This paper presents the results of the statistical analysis of biomass and NPP estimates from vast forest areas in relation to climatic conditions. The analysis was based on biomass data in 1236 sample plots as well as both NPP and biomass estimates in 480 sample plots that were established in mixed spruce and fir stands. Data were available for the entire Eurasian territory from the UK to South China and Japan, and included information on both total biomass and various biomass components. The climatic conditions of the sampled forests were characterized by using specific indices of zonality and continentality. A system of recursive equations was fitted to the data using standard regression analysis. For both the aboveground and total NPP, we found a statistically significant increase from the North to the South, but also a decrease from both the Atlantic and Pacific coasts towards the continental center of Siberia. We also found that lower upper-story NPP is compensated by higher understory NPP and vice versa. On the other hand, inconsistencies were identified concerning the share of root biomass in relation to the continentality index.


Introduction
Forests generally play an important role in the global carbon cycle (IPCC, 2013), which is an essential component of shaping the climate of the Earth (Schulze, 2000;Utkin, Zamolodchikov, & Milova, 2004;IPCC, 2013).The global carbon cycle is affected by net forest carbon balance, including deforestations causing massive emissions, and the net balance of existing forests and afforestations that altogether offset a significant portion of human-induced emissions.
Both the net carbon balance of the forests and the amount of carbon stored in the biomass and other forest carbon pools show large geographical and temporal variation (IPCC, 2013;Somogyi & Zamolodchikov, 2007) that are responses to variations in ecological conditions.Evidence is mounting that the net forest carbon balance can change in the short run (Fernando et al., 2014), has already changed due to climate change but also other factors (Spiecker, Mielikäinen, Köhl, & Skovsgaard, 1996;Somogyi, 2008), and may further change in the future especially due to climate change (e.g., Ludermilk et al., 2013).On large areas, forests may even become considerable sources of emissions in the future under certain climate change scenarios.In order to avoid or minimize adverse forest-related effects of climate change, it is imperative that we both have a better understanding of the effects of climate and climate change on the biomass and net primary production (NPP) of large forest areas, and mitigate climate change impacts by this understanding.
Both biomass and NPP are known to be linked to climatic patterns such as evapotranspiration (Rosenzweig, 1968;Lieth, 1974aLieth, , 1974b) ) as well as the amount of precipitation and average annual temperature (Lieth, 1974a(Lieth, , 1974b;;Luyssaert et al., 2007).Such links have often been studied indirectly by using geographical gradients (e.g.Shi, Sasa, & Koike, 2010).In this paper, we make an attempt to more directly relate forest biomass and productivity to climatic gradients.The analysis was done by using the Eurasian forest database, which was compiled and published by Usoltsev (2010).It includes forest biomass and NPP of sample plots together with consistent information on both climatic zonality and continentality of the geographical location of the stands.
Based on this unique set of data for the largest forest area in the world, i.e., the vast Eurasian forests, we found a statistically significant increase of both the aboveground and total NPP from the North to the South, but also a decrease from both the Atlantic and Pacific coasts towards the continental center of Siberia.We also found that lower upper-story NPP is coupled with higher understory NPP and vice versa.On the other hand, inconsistencies were identified concerning the share of root biomass in relation to the continentality index.

Data and Methods
We used all data of mixed spruce (Picea sp.) & fir (Abies sp.) forests from the database.Different species of these two genera usually grow together in the boreal zone, although with locally varying predominance of one or the other.A dedicated comparative study of spruce and fir-dominated forests in the Ural Mountains did not detect statistically significant differences in the biological productivity of these forest types between 20 and 130 years of age (Usoltsev, Vorobeichik, & Bergman, 2012), therefore, we used all data of such forests.The simultaneous analysis of the two different genera was also caused, in some cases, by the impossibility of separating them by eco-regions (e.g. when spruce and fir occur together in the boreal zone, or when Wilson's spruce and the Georgian fir occur together in the zone of mixed forests of China).In other cases, a combined analysis of different species was possible due to the fact that different but closely related species of the same genus, rather than the same species, occur throughout Eurasia (e.g.silver fir, Siberian fir or Chinese fir), demonstrating that the habitats of woody species are associated with specific eco-regions.The latter is known as vicariation in plant community chorology (Tolmachev, 1962) whereby ecologically equivalent species replace each other in the course of long-time geographical separation once they form a separate but continuous habitat.In order to explore the geographical distribution of woody NPP in the broadest geographical ranges, vicariation is inevitable to consider, and thus we included data from several vicariants in our database.
For this study, data of 1236 plots with estimates of biomass and data of 480 plots with estimates of both NPP and biomass were analyzed.The data that could be retrieved from the database and that included all relevant biomass fractions (stems, branches, foliage and roots of upper-story trees, and total biomass of the understory) are distributed in Eurasia as follows: in Europe, 744 plots in Norway spruce, silver fir and Siberian fir dominated stands; in Siberia and the Far East, 185 plots in Siberian spruce, Siberian fir, Jeddo spruce and Khingam fir dominated stands; in Northwest China, 77 plots in Asian spruce dominated stands; in Central and Southern China, 176 plots in Wilson's spruce, Georgian fir and Chinese fir dominated stands; and in Japan, plots in Jeddo spruce, Koyama spruce, Veitch's fir, Momi fir and Sakhalin fir dominated stands.Data on NPP and biomass in the sample plots are supplemented with mensuration characteristic of forest stands (such as tree age, density and growing stock).
We note that, in previous studies of the geography of birch, spruce-fir and Siberian (Korean) pine, data of only those plots were used where both NPP and biomass were estimated (Usoltsev, Noritsina, & Bornikov, 2010;Usoltsev, 2013).However, in order to obtain more reliable results regarding the geographical patterns over vast areas, we merged the two types of data sets in this study.
We suppose all plots were established, based on the Methodological guidelines for the International Biological Program (Minimum program, 1967), in typical "background" habitats, i.e. in representative locations of the spruce-fir plant communities.As the data were collected in many countries and under rather different ecological conditions over a longer period of time, it is assumed that the sample plots used are indeed representative and that thus they allow for an analysis of the geographic gradients of biological productivity.It is also assumed that all relevant biomass fractions were estimated following the above methodological guidelines.
The relevant biomass and NPP information from the selected plots was analyzed in relation to ecological variables.We selected two quantifiable variables to describe the relevant characteristics of Eurasia.One is based on a climatic classification, whereas the other is related to climatic continentality, which is estimated by using a specific type of continentality index.
One way to define climatic zonality is to use the sum of effective temperatures (SET) above + 5 °C or +10 °C as a continuous variable.However, the known scheme of SET by Tuhkanen (1984) is only established to the Northern part of Eurasia.Therefore, we used five broad climatic zones: (1) Subarctic, (2) Northern moderate, (3) Southern moderate, (4) Subtropical and (5) Subequatorial (Alisov & Poltaraus, 1974;Bazilevich & Rodin, 1967) where the boundary between Northern moderate and Southern moderate climatic zones was defined based on the Southern boundary of Southern taiga by Bazilevich and Rodin (1967).Figure 1 shows both the geographical distribution of these zones and the geographical locations of the sample plots within the zones.The above numbers are not continuous variables, however, such serial variables (e.g., Kraft's tree class from I to V, tree height class from I to X, site index from I to V etc.) were commonly used in modeling of forest dynamics (Tyabera, 1980;Yanovskii & Moiseev, 1985;Matveev & Usoltsev, 1991) as a semi-continuous proxy for modeling variables that could otherwise not be accounted for, and we treat them as such.
Figure 1.The location of the sample plots in the climatic zones: 1 -Subarctic, 2 -Northern moderate, 3 -Southern moderate, 4 -Subtropical and 5 -Subequatorial (Alisov & Poltaraus, 1974;Bazilevich & Rodin, 1967) There are a number of ways to quantify the degree of climatic continentality (Knoch & Schulze, 1952).Known formulas mainly differ in the relative contribution and/or ratios of two components: one, the amplitude of temperatures of the warmest and the coldest months, and two, the geographical latitude.In particular, the formula by Zenker (after Borisov, 1967) is as follows: where CI is the climatic continentality index, usually varying in the range from 8 to 100; А is the amplitude of annual temperatures, °C; and φ is geographic latitude, degrees.
In the study of the geographic patterns of larch biomass and NPP in Northern Eurasia, which analyzed 390 sample plots, the Zenker index of continentality (ranging from 35 to 95) was found to have significant negative correlation with the total and aboveground biomass and NPP (Usoltsev, 2001;Usoltsev & Koltunova, 2001;Usoltsev, Koltunova, Kajimoto, Osawa, & Koike, 2002).However, we used a somewhat more complex formula that uses the same variables and that was developed by Khromov (1957) based on the formula by Gorczynski (1920): Just like with the climatic zones, each sample plot was positioned according to the continentality index of their geographic location.This is demonstrated in Figure 2 using isoconts, i.e. lines of equal continentality index.
Figure 2. Allocation of the sample plots on the isocont map of Eurasia by Khromov (1957).Numbers on the map represent specific values of the Khromov continentality index (X axis is longitude, Y axis is latitude in degrees) An analysis based on the above types of variables inevitably concerns itself with phenomena of stochastic nature that involve correlation (or "co-relation" as proposed by F. Galton in 1888) between various variables.In order to identify how climatic conditions best explain patterns of biomass and NPP, standard multivariate regression analysis is used.However, regression was not applied to a single equation, rather, to a system of related ones.This was done by recursively fitting both the geographic/climatic information and biomass/NPP data to a system of related equations of the following general form: where N = tree density, 1000/hа; А = stand age, yrs; V = growing stock, m 3 /hа; Bi = biomass of i-th fraction (stems, Bs; branches, Bb; foliage, Bf; roots, Br; aboveground biomass, Ba; total biomass, Bt; and understory, Bu), t/hа; Zon = the serial number of climatic zones (Figure 1); CI = continentality index (Figure 2).
For the plots with NPP and biomass data, the following additional regression equation is applied: where Zi = annual NPP of i-th fraction: stem, Zs; branches, Zb; foliage, Zf; roots, Zr; aboveground biomass, Za; total biomass, Zt; and understory, Zu, t/hа.
Variables A, N and Zon have both a linear and a quadratic term in the equations, whereas only the linear term is used for all other variables.
Fitting Equations 3 and 4 was done in the following sequence: first, N was calculated from the values of A, Zon and CI; second, V was calculated using the values of N obtained in the first stage and the same values of A, Zon and CI; third, Bi was calculated using the values of N and V obtained in the first and the second stages and the same values of A, Zon and CI; finally, Zi was calculated using the values of N, V and Bi obtained in the three previous stages and the same values of A, Zon and CI.
When fitting the above equations, standard regression statistics were calculated.Note however that, due to the nature of the recursive system, the adjusted R 2 values reported always refer to one step of the entire procedure, i.e. to fitting an equation to values of its variables that are taken as independent ones in that step.
All variables in the system of Equations 3 and 4 were used after a logarithmic transformation.Such transformation was used because equations without log operators give negative values of biomass and NPP when the value of CI index is more than 80.

Results and Discussion
The regression analysis of the system of Equations 3 and 4 produced both the required regression coefficients and the coefficients of determination, adjR 2 (Tables 1 and 2).All regression coefficients of the independent variables are significant at the 95% probability level.In order to demonstrate the relationship between NPP and both the climatic zones and the continentality index, we calculated the values of NPP fractions from the equations at specific reference values of age, climatic zone and CI. Figure 3 shows NPP values at the age of 100 years and for climate continentality index of 75 by climatic zone, demonstrating low NPPs towards the North and higher NPPs towards the South.Figure 4 shows the effect of continentality at a fixed age of 100 years and in the climatic zone: all NPPs decrease as continentality increases.Finally, Figure 5 demonstrates that NPP of the understory, estimated for similar specific reference points using the above procedure, also show strong but sometimes different relationships with climate and continentality: with the exception of the northernmost areas, and just like upper-story NPP, understory NPP generally increases towards the South, whereas it increases with increasing continentality, which is the opposite for upper-story NPP.
Upper-story NPP, t•ha -1 •yr -1 Climatic zone Consistent with the above, Shi, Sasa, and Koike (2010) found negative trend of 30-year-old larch forest biomass over latitude (the latitude series of 47, 50, 52 and 62° N correspond to Southern, Central and Northern part of the Daxingan mountains as well as of Central Siberia).According to their findings, changes in latitudes from 47 up to 62° N are correlated with the drop of forest biomass to 1/3 -1/4 of its southernmost value.It is also shown that forest NPP is decreasing synchronously in the direction from the tropics to the two poles (Anderson, Allen, Gillooly, & Brown, 2006).For both forest biomass and NPP, a statistically significant positive relationship was found for spruce, birch, natural Scots pine stands and Scots pine plantations in Northern Eurasia over the average monthly sum of effective temperatures above + 5 °C (an index developed by Tuhkanen, 1984), within the range of the index between 20 to 80 °C, (Usoltsev, 2007;Usoltsev et al., 2010).The dependence of forest productivity on latitude (i.e., the location along a North-South direction) is thus reconfirmed.
It should be noted that the isothermal lines, which can also be taken as the boundaries of natural climatic zones, are correlated with the latitudinal calibrating of the Earth's surface, although with some deviations due to the altitudinal zonality of mountain territories.Therefore, the latitudinal gradient of forest biomass and NPP can be described quantitatively in relation either to the sum of effective temperatures or the serial number of a natural zone in the direction from North to South or from South to North, and directly in relation to latitude (Lavrenko, Andreev, & Leont'ev, 1955;Usoltsev, 1998).
It is more complex to relate forest productivity to the meridional gradient.Komarov (1921) developed the concept of meridional zonation of plant cover that complements the latitudinal zonation and should be taken into account in the allocation of biogeographical regions.Komarov (1921) also distinguished two types of vegetation on the major continents: the seaside type, situated in a narrow strip along the coasts, and the inland one, away from the seaside type.When intersecting with seven latitudinal zones, they form 42 broad floristic districts, each with its own climate, soil, their plant endemism and the predominant type of vegetation cover.Thus, the main plant cover variability occurs not only latitudinally due to a change in the intensity of solar radiation from South to North, but also meridionally (from the sea coasts to the center of a continent) as a result of changes in continentality (i.e., changes of both temperature and moisture conditions, Volobuyev, 1947;Kurnayev, 1973;Nazimova, 1995).Summarizing the extensive more recent literature, findings with respect to forest productivity going from the West to the East of Russia are contradictory: according to Pobedinskii (1964), it decreases, and according to Tyabera (1988), it increases.Our study, however, has not produced estimates along a simple geographical direction, rather, along an ecological gradient, i.e. continentality, and found a rather strong dependence of forest productivity on continentality.
The NPP estimates of spruce-fir forests are significant over both climatic gradients, however, some NPP fraction demonstrate different patterns.Stem NPP as well as aboveground and total NPP increase monotonically from the 1st to the 5th climatic zone, with a ratio of the NPP in the climatic zone 5 relative to the NPP in climatic zone 1 being 11.5 and 6, respectively.Branch, foliage and root NPP over climatic zones show a bell-shaped curve, reaching a maximum in climatic zones 3 and 4 (Figure 3).
In climatic zone 3, the NPP of all fractions decreases monotonically from the Atlantic and Pacific coasts towards the center of continentality but the rates of this decline are rather different: the branches show the sharpest decline (the NPP of the branches in areas of the highest CI is only 1/28th of that in areas of the lowest CI), whereas the foliage shows the smallest difference of 30% of NPP from areas of low CI to areas of high CI (Figure 4).This is consistent with earlier findings: under self-thinning conditions, the biomass of basic forest-forming species in Northern Eurasia showed a statistically significant decline in the direction from the South to the North (related to the decreasing SET), and from the Atlantic and Pacific coasts to the continentality center in Siberia (Usoltsev, 2003).The biological and ecological causes of the latter are not yet fully understood.
Going from the South to the North, the share of roots in the total biomass (see also Usoltsev, 2001Usoltsev, , 2003) ) and in the total NPP (Figure 3) increases from 22 to 42% and from 23 to 35%, respectively.This increasealso seems to be related to the decline of SET, although other factors such as soil properties might also be important.In contrast, we have found an inconsistent pattern concerning roots over continentality: the share of the roots in the total biomass increases from 23 to 33% (Usoltsev, 2003) in direction from the Atlantic and Pacific coasts to the continentality center in Siberia, but the share of the roots in the total NPP decreases from 35 to 17% in the same direction (Figure 4).Currently, it is not possible to explain the causes of this inconsistency.
Concerning the methodology applied, it must be considered that forests are complex biological dynamic systems that often require multidimensional analysis of its structure and dynamics.When modeling basic processes in forest ecosystems, using a number of partly or completely interdependent relationships may often produce results of acceptable accuracy.Mathematically, such relationships may take the form of a system of related equations that can efficiently deal with inherent variability.We used such a simultaneous system of regression equations so that we fitted one equation at a time that was followed by recursively fitting others in a logical sequential order (Amateis, Burkhart, Greber, & Watson, 1984).
Currently, many phenomena and features in forests can hardly be estimated, let alone measured, and this has a bearing on the mathematical models used.However, this limitation is technical rather than one of principle, and its elimination is a continuous process.The method of factorial design of experiments (e.g., the rather standard method of passive factorial design by Nalimov, 1971) suggests the use of not only quantitative but also of qualitative factors which may be quantified and used as algorithms of regression analysis (Kleijnen, 1978).
In general, it is assumed in regression analyses that the variables are continuous.However, "if the factors are qualitative, then purely mnemonic numbers are attached to their levels" (Kleijnen, 1978, p. 21)."For example, colors or even shade may be expressed in numbers only because we know the rainbow" (Ivakhnenko, 1969).Adler, Markova, and Granovskii (1976) suggest that "the boundary between the concepts of a qualitative and a quantitative factor is very nominal".Aivazyan, Enyukov, and Meshalkin (1985) considered several types of variables: quantitative, ordinal (numerical and ordered), nominal (numbered but not ordered), and recommended appropriate mathematical tools.Kleijnen (1978), on the other hand, who paid much attention to the problem of the robustness of the evaluation in multi-factor experiments, does not see the need to introduce any limitations (contraindications), at least for the first two types of variables.
Admittedly, when discontinuous variables are used, there is "a model with errors in the variables" (Seber, 1980).In these cases, a specific computation algorithm is required, leading to a more complex method than the usual least-squares method.However, according to Seber (1980), if the values of the variables may be tested, the model can be applied as if the variables were devoid of errors.
Another limitation of multiple regression models is the multicollinearity (i.e., mutual correlation) of the factors that must be orthogonal, i.e. they must act independently of each other.The two geographic/climatic factors that we used in our analysis, i.e. the natural zonality and continentality index, designating correspondingly SET and the amplitude of annual temperatures, are orthogonal in principle, at least if we use the continentality index by Polozova (1954).This index is to be calculated for any given latitude as it does not involve the factor of geographic latitude by using the following form: where А pos = maximum positive anomaly of air temperature, °C; A neg = maximum negative anomaly of air temperature, °C; a i = anomaly of air temperature, °C.
However, as it was shown by Usoltsev (2003), the index by Polozova is less acceptable than the index by Zenker because it shows continentality index equal to 70 both for Siberia and for Japan that in fact have rather different continentality.The factors used in our study are partly orthogonal because the index by Khromov uses the amplitude of temperatures and the geographic latitude.The latter is correlated with SET.The formula by Khromov and the corresponding scheme of isoconts were the most preferable because the schemes of isoconts fitted according to formulas by Zenker (after Borisov, 1967) and by Conrad (1946) were only extended on the territory of Northern Eurasia (Borisov, 1967;Tuhkanen, 1984).
A system of recursive equations needs to apply some specific calculation methods, i.e. two-stage, three-stage and multistage least squares (Furnival & Wilson, 1971;Borders & Bailey, 1986;LeMay, 1990) to exclude correlations of both coefficients and residuals in the subsequent fittings.However, this was not tested in Equations 3 and 4.But we assume that these correlations are low at a high probability as it was reported by Amateis et al. (1984) when analyzing the recursive system of three equations with 5 exogenous variables and 3 endogenous ones: "The recursive system of equations was refitted to the combined data and the variance-covariance matrix of residuals was tested to determine if the off-diagonal elements were significantly different from zero using a test adapted from Anderson (1958).The results showed that the assumption of independence between equations in the recursive system is not unreasonable." It is to be further considered that any practical benefits in ensuring the correctness of a mathematical model is primarily achieved by an appropriate database and a correct underlying theory rather than a statistically sound procedure of parameter estimation.In practice, the logical validity of the model matters more than accuracy.The choice of more powerful means of mathematical statistics may not be a panacea (Klein, 1960).For example, using harvest data of root biomass-height profiles in Pinus sp., the usual multi-factorial regression model, characterized by an R 2 of 0.97 and by estimation error 9.8% was designed (Usoltsev & Krepkii, 1994).Upgrading the model by preserving statistical correctness (i.e., by using normality test by Ghosh (1996) and bias correction following Finney (1941), and Akaike's (1974) criterion), the same root biomass data gave an R 2 of 0.40, and an estimation error of 27% (Hoffmann & Usoltsev, 2001).Thus, when designing regression models, their statistical correctness may not suffice, rather, model structure correctness should be ensured by means of preliminary analysis of the initial data (Liepa, 1980).The theory behind our approach strongly support the model structure selected.This, and the high quality of the data in the database also support the high validity of the conclusions.
The results of our analysis can be useful, among others, in calculating the input part of carbon balance in forests that can significantly contribute to mitigating climate change.The equations reported and the regression statistics can help validate various process models and phenomenological simulations (Wolf, Ciais, Bellassen, Field, & Berry, 2011), and they can be used in forest resource assessments and forest economics as well as in forest carbon projections.However, further research might be needed e.g.concerning how the established relationships could change as a result of climate change.

Figure 3 .NPPFigure 4 .Figure 5 .
Figure 3. Estimates of the NPP of stems (a), branches (b), foliage (c), roots (d), aboveground (e) and total biomass (f) over climatic zone in spruce-fir stands of 100 years of age and of continentality index of 75

Table 1 .
Regression results of fitting the system of Equations 3 * Number of data points.

Table 2 .
Regression results of fitting Equations 4