Carbon Balance and Management

BioMed Central

Research

Open Access

Simulating the carbon balance of a temperate larch forest undervarious meteorological conditions
Motomu Toda*†1, Masayuki Yokozawa†2, Akihiro Sumida1,
Tsutomu Watanabe3 and Toshihiko Hara†1
Address: 1Biosphere Dynamics Research Group, Institute of Low Temperature Science, Hokkaido University, Sapporo 060-0819, Japan,
2Department of Global Resources, National Institute for Agro-Environmental Sciences, Ibaraki 305-8604, Japan and 3Cryosphere Environment
Research Group, Institute of Low Temperature Science, Hokkaido University, Saspporo, 060-0819, Japan
Email: Motomu Toda* - todam@pop.lowtem.hokudai.ac.jp; Masayuki Yokozawa - myokoz@niaes.affrc.go.jp;
Akihiro Sumida - asumida@hassaku.lowtem.hokudai.ac.jp; Tsutomu Watanabe - t-wata@lowtem.hokudai.ac.jp; Toshihiko Hara - thara@pop.lowtem.hokudai.ac.jp* Corresponding author †Equal contributors

Published: 30 May 2007
Carbon Balance and Management 2007, 2:6

doi:10.1186/1750-0680-2-6

Received: 24 January 2007
Accepted: 30 May 2007

This article is available from: http://www.cbmjournal.com/content/2/1/6© 2007 Toda et al; licensee BioMed Central Ltd.
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0),
which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract
Background: Changes in the timing of phenological events may cause the annual carbon budgetof deciduous forests to change. Therefore, one should take such events into account whenevaluating the effects of global warming on deciduous forests. In this article, we report on the
results of numerical experiments done with a model that includes a phenological module simulating
the timing of bud burst and other phenological events and estimating maximum leaf area index.
Results: This study suggests that the negative effects of warming on tree productivity (net primary
production) outweigh the positive effects of a prolonged growing season. An increase in air
temperature by 3°C (5°C) reduces cumulative net primary production by 21.3% (34.2%). Similarly,
cumulative net ecosystem production (the difference between cumulative net primary production
and heterotrophic respiration) decreases by 43.5% (64.5%) when temperatures are increased by
3°C (5°C). However, the positive effects of CO2 enrichment (2 × CO2) outweigh the negative
effects of warming (<5°C).
Conclusion: Although the model was calibrated and validated for a specific forest ecosystem, theimplications of the study may be extrapolated to deciduous forests in cool-temperate zones. These
forests share common features, and it can be conjectured that carbon stocks would increase insuch forests in the face of doubled CO2 and increased temperatures as long as the increase intemperature does not exceed 5°C.

BackgroundNowadays, there is much apprehension about the great
effects of global warming on the forest ecosystems, especially at mid to higher latitudes [1]. The abrupt climatechange with the effect of global warming may change distributional pattern of plant species through competition
among plants, and the change may then affect the global

carbon cycle [2] because the vegetation ecosystems, especially forest ecosystems play a significant role in carbon
exchanges between an atmosphere and biosphere [e.g.,[36]].
Global warming may also affect phenology of trees, especially deciduous trees. Then, the time changes in phenoPage 1 of 18(page number not for citation purposes)

Carbon Balance and Management 2007, 2:6

logical events may bring about changes in the annual
carbon budget in forest dynamics of the deciduous trees[7-9]. Phenology reflects the eco-physiological responsesof trees that are changed by daily forcing factors such as air
temperature, solar radiation and drought conditions in
the atmosphere and soil. Major phenological events ofdeciduous woody species are the timing of bud burst (leafunfolding), flowering and leaf-fall through leaf senescence with the course of seasons, i.e., the onset and end,
respectively, of the photosynthetic active period [10-13].
A change in phenological events might be one of the mostimportant factors affecting the future carbon sequestration potential of forests ([14,15]). In terrestrial land surfaces, percentages of deciduousness show a unimodal
distribution pattern with a peak at mid-latitude, while abimodal distributional pattern with two peaks, one atlower and the other at higher latitudes, is observed for thepercentages of evergreeness [16]. However, deciduoustrees are widely distributed over the temperate to boreal
climate zones, adapting to the local environmental conditions (e.g., Larix dahurica forest stand in eastern part ofSiberia). Thus, to elucidate the effects of global warmingon the deciduous forests in these climate zones and therole of deciduous forests in the global carbon cycle, it isrequired to develop a phenological model for deciduoustrees [17].
Meanwhile, the annual leaf area is also one of the mostsignificant determinants to carbon budget. Thus we havean interest in what determines the leaf amount each year,
especially maximum leaf area. The determination of maximum leaf area is closely related to phenological processes, based on the carbon reserves, and should be,therefore, considered in conjunction with the phenological events. However, an attempt has rarely been made tomodel the annual production of leaves (i.e. leaf areaindex, LAI) in deciduous trees (but see [18,19]). In mostof previous atmosphere-biosphere interaction models,
LAI has been given as one of assigned constants or external forcing data in relation to plant functional types (PFT)(e.g. [20-22]). For instance, Foley et al. [23] obtained LAI
of each PFT by simply dividing leaf carbon by specific leafarea, and Sitch et al. [2] represented LAI by dividing leaf
mass by crown area which was obtained from the allometric relationship between crown area and stem diameter.
Sato, Ito and Kohyama [24] determined the maximum
leaf biomass for each individual with three constraints
such as crown surface area, cross-sectional area of sapwood and available resources. In addition, Ito andOikawa [25] employed time-varying LAI so as to produce
maximum daily net carbon uptake from the daily grossprimary production at the entire forest stand level, whichis the optimal LAI concept proposed by [26]. However,such treatments of LAI do not consider trees' physiologicalresponses to varying environmental conditions. Aber et al.

http://www.cbmjournal.com/content/2/1/6[18] developed a model which estimates LAI in a year
using carbon reserves accumulated in buds and woody
organs during several years before the corresponding year,
and incorporated the framework into a forest carbon andwater balance model (PnET). It might be important toconsider seasonal carbon reserve mobilization for determining LAI of deciduous trees from an ecological point ofview because the dynamics of carbohydrate in tree species
are affected by annual variations in climate [27].
We developed a deciduous canopy development model torepresent seasonal and annual dynamics of deciduoustrees. The canopy development model includes a phenological model to determine the timing of bud burst andthe end of the leaf appearance and a model for estimating
maximum LAI. In this paper, we examined the population
dynamics and carbon balance of deciduous trees by coupling the canopy development model with an integrated
climate and forest dynamics model, the Multilayered Integrated Numerical Model of Surface Physics-Growing
Plants Interaction (MINoSGI) [28]. MINoSGI is one ofintegrated atmosphere-forest dynamics interaction models. The model describes population dynamics of trees or
growth (carbon sink) and death (carbon source) of individual trees which are relevant to carbon dynamics in aforest by using a partial differential equation for the size
distribution density [28-30]. Thus, MINoSGI is an alternative approach to individually based vegetation dynamics
models which resolve spatial locations of individual
plants to mimic horizontal interference between individuals [e.g., [31,24]]. Here, the original algorithm of
MINoSGI, which had been applicable only to evergreen
trees, was further extended to be applicable to deciduoustrees by incorporating the above-mentioned features withrespect to the timing of bud burst and the seasonal carbon
storage dynamics of deciduous trees. The model wasapplied to a deciduous coniferous Larix kaempferi forest ina climatologically temperate area, and its performancetested by comparison with inventory data.
Further, sensitivity analysis was conducted to assess anindividual or combined effect of increased CO2 concentration and air temperature on forest dynamics and production of a forest over a hundred years. The objective is to
examine the responses of size-structure dynamics of aplant community to changes in exogenous factors such as
CO2 concentration, air temperature. We emphasize that it
is not intended to project impacts of future climate change
on plant communities and global carbon cycle. Therefore,
the sensitivity analysis is the first step to evaluate theimpact of predicted changes in climate conditions such as
air temperature and atmospheric CO2 concentration withthe effect of global warming on forest production alongwith variations in size structure of the deciduous coniferous forest on a long time scale.

Page 2 of 18(page number not for citation purposes)

Carbon Balance and Management 2007, 2:6

http://www.cbmjournal.com/content/2/1/6

Simulated biomass, averaged tree height and organic
matter production
It was examined whether the model with the adjusted
parameters (see the Methods section in details) produces
realistic population dynamics of the forest over sevenyears. The simulated annual variations in the aboveground biomass and averaged tree height were reproduced quite well by simulation (Fig. 1). The observed andsimulated biomass and averaged tree height were 113,
111 ton ha-1 and 14.2, 14.5 m in 1989, respectively.

(a)

(a)

8

6

0.02

4
0.01
2
0
0

120

1982 1983 1984 1985 1986 1987 1988 1989 1990

100
80
60

14

(b)

Rleaf (kg m-2 day-1)

0.02

1982 1983 1984 1985 1986 1987 1988 1989 1990

Averaged tree height,
have (m)0.03

LAI observed
LAI simulated
GPP simulated
NPP simulated

Air temperature
Rleaf simulated

(b)

0.015

290

0.01
280
0.005
270

0
-0.00512300

Air temperature (K)Biomass (ton ha-1)

140

Observed
Simulated

Sensitivity to meteorological conditions
Changes in the future climate system might affect the
dynamics of forests that are generated by competition forresources between individual trees. Then, dynamicchanges of forests have a possibility to produce variations
in carbon balance of the forests. Herein, with the fully validated MINoSGI model, we aimed at sensitivity analysis
to assess how changes in air temperature and atmospheric

LAI (m2 m-2)

The simulated annual gross primary production (GPP)and net primary production (NPP) within the forest in
our study over 7 years were 32.6 and 12.1 on average inthe range from 31.8 to 33.9 and from 10.8 to 14.7 (ton ha1 year-1), respectively (Fig. 2(a)). Both GPP and NPP havea similar seasonal trend, i.e., two peaks in spring andautumn. A large decrease in GPP and NPP was shown insummer. The decrease was caused by the increase in leaf
respiration (Rleaf) which is mainly controlled by a seasonal change in air temperature (Fig. 2(b)). Rleaf was seventimes as large as woody organ's respiration rates (stem androot) and 1.8 times larger than soil respiration on average.
Therefore, averaged NPP/GPP was 0.37. In addition, the

NPP simulated in the present study was consistent withthe results referred to in [32], who summarized NPPs estimated in many Larix forests in the cool-temperate zone inChina, and was equivalent to those obtained in otherdeciduous forests (e.g., 10 ton ha-1 year-1 in the case of amixed forest which consisted of sugar maple, American
beech, and yellow birch in [33]).

GPP, NPP (kg m-2 day-1)

Results

260
1982 1983 1984 1985 1986 1987 1988 1989 1990

Time (year)
1081982 1983 1984 1985 1986 1987 1988 1989 1990

Time (year)
biomass
Observedkaempferi
Figure
1and
stand
and
(b)simulated
averagedannual
tree heightvariations(have)inin(a)theforest
Larix standObserved and simulated annual variations in (a) forest standbiomass and (b) averaged tree height (have) in the Larix
kaempferi stand.

Figure(a)production
intheSimulated
Larix
2 kaempferi(GPP)
seasonal
andstandnet
andprimaryannual variations
productionin (NPP)gross primaryand LAI(a) Simulated seasonal and annual variations in gross primary
production (GPP) and net primary production (NPP) and LAI
in the Larix kaempferi stand. The vertical axes on the left- andright-hand sides are employed for GPP, NPP and LAI, respectively. LAIs observed in 1982, 1983 and 1984 are also plottedin this figure. (b) Simulated seasonal and annual variations inleaf respiration (Rleaf) and air temperature. The vertical axeson the left- and right-hand sides are employed for Rleaf and airtemperature, respectively.

Page 3 of 18(page number not for citation purposes)

Carbon Balance and Management 2007, 2:6

CO2 concentration affect production and size structure ofa L. kaempferi forest with forest development. In the simulation, we assumed a monospecific L. kaempferi forest
which was planted at a density of 5000 trees per hectre.
The initial sapling size in the stand was assumed to be uniform and was set at the smallest height class at the beginning of the simulation. In addition, we used an annual
meteorological dataset with 150 repetitions, which was
generated using the data in 1986 from the Inabu experimental site ([34]; Fig. 3) in order to exclude the effect ofyear-to-year variations in meteorological forcing data on
production of the Larix forest. In addition, we assumed
that the values of initial carbon stock in each component
in the soil carbon dynamics model were identical to thoseused in the model validation for the Larix forest site inInabu experimental field [See the Methods section]. Eightsimulations were made under the following conditions:(1) the present climate conditions (control; see Fig. 3), (2)3°C increase in air temperature (+3°C), (3) 5°C increasein air temperature (+5°C), (4) 10°C increase in air tem-

-1) data
-2),
(K),
the
FigureAMEDAS
specific
3variations
humidity
meteorological
(kgkg
and
precipitation
station(b)ofairInabu,temperature
obtained
Seasonalmeanbasis;(a)solarinradiation
meteorological(WmvariablesonaJapandaily atSeasonal variations in meteorological variables on a dailymean basis; (a) solar radiation (Wm-2), (b) air temperature(K), specific humidity (kgkg-1) and precipitation obtained at
the AMEDAS meteorological data station of Inabu, Japan. The
annual dataset with 150 repetitions was used as part of meteorological forcing data.

http://www.cbmjournal.com/content/2/1/6

perature (+10°C), (5) doubled CO2 concentration (2 ×CO2), (6) combined effect of (2) and (5), (7) combined
effect of (3) and (5), (8) combined effect of (4) and (5).
The annual averaged values of daily mean solar radiation,
air temperature, specific humidity and annual amount ofprecipitation were 145 Wm-2, 7.7°C, 5.4 kgkg-1 and 2256
mm, respectively (Fig. 3). The precipitation accumulated
during summer (from July to September) was 1196 mm,
representing more than a half of the annual amount (Fig.
3).
We can see similar annual trends of slightly increasing or
almost constant GPP and decreasing NPP and NEP (netecosystem production) with the decrease in tree density,
namely, with forest development, except for the one with
the individual effect of increased air temperature by 10°C
(Fig. 4(a)–(c)). In particular, a large decrease in annual
NPP with forest stand development occurred in the sevensimulations due to a large increase in autotrophic respiration rate (Ra) relative to GPP with tree growth in 2 × CO2
condition. Ra was increased with forest development in
the seven simulations due to the increment of sapwood
and root biomass with tree growth from seedling (Fig.
4(d)). In contrast, heterotrophic respiration rate (Rh) with
forest development got greater as air temperature
increased (Fig. 4(e)). Rh is more strongly influenced by thelitter decomposition pool with greater decomposition
rate than by the remainder soil organic matter (SOM)
pools with smaller decomposition rate. The amount of litter supply from above-ground into the litter decomposition pool appears to be closely related to the size of trees.
LAI gradually decreased with forest development in all theeight simulations. LAI at a specific tree density (i.e. forest
development) decreased in the order of 2 × CO2, +3°C
and 2 × CO2, +5°C and 2 × CO2, control, +10°C and 2 ×CO2, +3°C and +5°C, the same as in GPP, NPP and NEP.
The great decrease in LAI under lower tree densities wascaused by the individual effect of increased temperature
which yielded a number of dead trees (Fig. 4(f)).
GPP, NPP and NEP in 2 × CO2 were greater than those inother conditions. 2 × CO2 enhanced the growth of trees
with forest development and produced the greatest averaged tree height, stem diameter and LAI in all the conditions (Fig. 4(f)) and then brought about increases in the
cumulative NPP and NEP by 55.9% and 53.2%, respectively, over the 150-year period, compared with thoseunder the control conditions. Furthermore, the decreases
in annual NPP and NEP with forest development whentree density was less than 0.1 were the largest in 2 × CO2,
and annual NPP and NEP after 150 years from the start ofthe simulation were 1.24 and 0.15 kgm-2. These values
represented decreases by 30% and 74% as compared withthose obtained at the final time step when tree density was
assumed to be equal to the initial one, which is identical

Page 4 of 18(page number not for citation purposes)

Carbon Balance and Management 2007, 2:6

(a)
4

(c)

(b)

2
1.5

2

(8) +10oC, 2xCO2

NPP (kgm-2)

GPP (kgm-2)

3.5

2.5

(5) 2xCO2(6) +3oC , 2xCO2(7) +5oC , 2xCO2(1) control(2) +3oC(3) +5oC(4)+10oC

3
2.5
2
1.5
1

NEP (kgm-2)

4.5

http://www.cbmjournal.com/content/2/1/6

1.5
1

1
0.5

0.5

0

0

-0.5

0.5
0
0.5

0.4

0.3

0.2

0.1

0.5

0

0.4

0.3

0.2

0.1

0

0.5

0.4

0.3

0.2

0.1

0

0.1

0

Tree density (m-2)
0.5

Heterotrophic
respiration rate (kgm-2)

Autotrophic
respiration rate (kgm-2)

1.6
1.2
0.8
0.4
0
0.5

0.4

0.3

0.2

Tree density (m-2)

0.1

0

(f)

(e)

10

0.4

8

LAI (m2m-2)

(d)
2

0.3
0.2

6
4

0.1

2

0

0
0.5

0.4

0.3

0.2

Tree density (m-2)

0.1

0

0.5

0.4

0.3

0.2

Tree density (m-2)

index
Figure(LAI)
4 annual
with thedecreasein treedensity (forest
standprimarydevelopment)(e) heterotrophic
respiration
rateproduction(Rh) and (f)(NPP),
leaf area
ecosystem
Simulated
production
variations(NEP),
in forest(d)autotrophic
production;
respiration(a) gross
rate(Ra),production(GPP), (b)net primary(c) netSimulated annual variations in forest production; (a) gross primary production (GPP), (b) net primary production (NPP), (c) net
ecosystem production (NEP), (d) autotrophic respiration rate (Ra), (e) heterotrophic respiration rate (Rh) and (f) leaf areaindex (LAI) with the decrease in tree density (forest stand development). Eight simulations were conducted under the following conditions: (1) the present climate conditions (control; see Fig. 7), (2) 3°C increase in air temperature (+3°C), (3) 5°C
increase in air temperature (+5°C), (4) 10°C increase in air temperature (+10°C), (5) doubled CO2 concentration (2 × CO2),(6) combined effect of (2) and (5), (7) combined effect of (3) and (5), (8) combined effect of (4) and (5).

to the output from static land surface models that do notconsider forest dynamics.
In 2 × CO2, the growth of all trees was enhanced and alarge difference in the growth rate between large- andsmall-sized trees was also induced (Fig. 5). During theperiod with tree densities greater than 0.1, intensive treecompetition reduced tree density with time, and a numberof smaller-sized trees in the stand died gradually due tothe light shade effect of larger-sized ones. However, at tree
densities less than 0.1, the relative growth rate of largersized trees in the forest was greatly reduced by thedecreased carbon assimilation essential to their growth
and the increased Ra (Fig. 4(d)). Further, Rh was decreased
during the corresponding period due to a slight supply ofleaf-fall from living larger-sized trees. Hence, we see
deaths of larger-sized trees and the great reduction ofannual NPP and NEP. This was because of the increase invertically overlapped leaves between smaller- and larger-

sized trees as a consequence of large height-growth ofsmaller-sized trees. During the later simulation period,
the height of larger-sized trees in the stand exceeded 50 m.
Accordingly, due to the deaths of larger-sized trees the
reduction of annual NPP and NEP in 2 × CO2 was thegreatest in all the conditions.
On the other hand, an individual effect of increased airtemperature (+3°C, +5°C, +10°C) brought about
decreased NPP and NEP compared with the control conditions. The increase in air temperature brought about
earlier start of bud burst and thus the increase in the
period of photosynthetically active leaves. In +3°C and
5°C bud burst was earlier by 15.6 days (97.7 days from
Jan. 1) and 26.8 days (86.4 days from Jan. 1) than that in
the control conditions (113 days from Jan. 1) on averagethroughout the simulation period. Nevertheless,increased temperature did not enhance the growth oftrees, especially the growth rate of trees in +10°C was

Page 5 of 18(page number not for citation purposes)

Carbon Balance and Management 2007, 2:6

0.3

http://www.cbmjournal.com/content/2/1/6(a) U=0.5

0.3(5) 2xCO2(1) control(2) +3oC(3) +5oC(4) +10oC(b) U=0.4

0.2(c) U=0.3(6) +3oC , 2xCO2(7) +5oC , 2xCO2

0.16(8) +10oC, 2xCO2

0.2

0.2

0.1

0.1

0.12
0.08
0.04

0

0
0

5

10

15

20

Frequency (m-2)(d) U=0.2

0.1

0.12

0
0

5

10

15

20

(e) U=0.1

0.03

5

10

15

20

25

40

45

(f) U=0.05

0.08
0.02

0.06

0.08

0.04
0.04

0.01

0.02

0

0
5

0.02

0

10

15

20

25

(g) U=0.04

30

0.016

0.016

0
10

15

20

25

30

35

(h) U=0.03

15

0.01

20

25

30

35

(i) U=0.02

0.008

0.012

0.012

0.006
0.008

0.008

0.004
0.004

0.004
0

0.002

0
15

20

25

30

35

40

45

50

0
15

20

25

30

35

40

45

50

15 20 25 30 35 40 45 50 55 60

Tree height, h (m)FigureSizeρ
= 0.02
distributions
5
obtained from the eight simulations in Fig. 4 at 9 specific tree densities from the initial tree density ρ = 0.5 toSize distributions obtained from the eight simulations in Fig. 4 at 9 specific tree densities from the initial tree density ρ = 0.5 toρ = 0.02.
much smaller than that in other conditions (Fig. 5) and all
the trees died at year 123 from the start of the simulation
period. This suggests that the negative effect of an increase
in respiration by increased temperature on NPP and NEP
dominates over the positive effect of prolonged growing
season on NPP and NEP. Hence, the magnitude andannual variations in NPP and NEP were reduced byincreased temperature (Fig. 4(b), (c)). On the other hand,
Rh decreased as air temperature increased, and Rh did notincrease with forest stand development because trees withsmaller growth rates did not much contribute to the litter
decomposition pool. The increase in air temperature by
3°C or 5°C brought about the reductions of 21.3% and34.2%, 43.5% and 64.5% in cumulative NPP and NEP,
respectively, over the 150-year simulation period, compared with those under the control conditions. Furthermore, cumulative NPP decreased by 440% in +10°C
compared with the control conditions, and cumulative
NEP was negative, indicating that the forest worked as carbon source.

The combined effects of increased temperatures and doubled CO2 concentration produced intermediate production between those in each individual effect, and NPP got
smaller as air temperature increased in 2 × CO2 conditions. In 2 × CO2 and +3°C, cumulative NPP and NEP
increased by 35.2% and 11.9% over the 150-year simulation period, compared with those under the control conditions. These increases in NPP and NEP indicate that thepositive fertilization effect of 2 × CO2 on NPP and NEP
dominated over the negative effect (due to increased respiration) of increased air temperature by <5°C. In contrast, cumulative NPP increased by 18%, but cumulative
NEP decreased 14.4% in 2 × CO2 and +5°C comparedwith the control conditions. Furthermore, cumulative
NPP and NEP decreased by 33.4% and 109.5% in 2 × CO2
and +10°C compared with the control conditions. This
indicates that the negative effect of increased air temperature by ≥ 5°C on NEP dominated over the positive effectof 2 × CO2.

Page 6 of 18(page number not for citation purposes)

Carbon Balance and Management 2007, 2:6

Discussion
A lack of the understanding of plant acclimation to futureclimate change makes it difficult to elucidate the role ofthe forest in global carbon cycle on a long time scale. Furthermore, the lack produces uncertainties of the sensitivity tests of plant responses to doubled CO2 concentrationand increased air temperature conditions. Recently,
Körner [35] has comprehensively summarized the resultsof plant CO2 responses obtained from CO2-enrichment
experiments and represented that elevated CO2 concentration is beneficial for plant growth irrespective of species
and experimental designs. Some experimental reports
showed that elevated CO2 concentration stimulated
strongly only initial growth of trees only, whereas suchplant CO2 responses declined or disappeared after several
years [36,37]. The difference in plant CO2 responses withage has probably arisen from acclimation of trees to theelevated CO2 concentration. Our sensitivity tests did notconsider the plant acclimation. Therefore, the predicted
results may be different from the above-mentioned experimental ones concerning the impacts of elevated CO2 concentration and increased air temperature on the
production of a deciduous forest.

Conclusion
It is well-known that rising CO2 concentration has beencaused by anthropogenic emissions. Simultaneously, theclimate change due to the rising CO2 concentration
should accelerate carbon emission from natural forest
ecosystems into the atmosphere [38] which means positive feedback of the release of carbon into the atmosphere.
It was suggested that global carbon cycle is strongly controlled by soil decomposition processes of soil organiccarbon and it is important to fully consider soils in assessing the full impact of climate change through a terrestrial
carbon cycle modeling approach [39,40]. Accordingly,
forest modeling is expected to be used as a tool for policy
decision-makers in order to predict the effect of climatechange on carbon stocks in the forest ecosystems and toplan the most sustainable management under climatechange [41].
For this purpose, a canopy development model for phenological processes of the timing of bud burst and for theestimation of annual variations in leaf amount (LAI) by
considering seasonal carbon reserve was designed. Thismodel was coupled with the integrated climate and forest
dynamics model (MINoSGI), and its application to a
deciduous L. kaempferi forest was examined for stand-levelvariables such as biomass, averaged tree height and NPP.
Although the present model deals with averaged sizedependence of individual growth in a limited one-dimensional framework, our model simulation yielded realistic
forest dynamics over seven years as the result of carefulparameter adjustment based on the inventory data.

http://www.cbmjournal.com/content/2/1/6

The agreement between observed and simulated values ofbiomass and averaged tree height enhanced the reliabilityof MINoSGI, whereas a discrepancy between the observed
and simulated timing of bud burst in 1983 occurred. Several species-specific parameters relevant to the timing of
bud burst in the phenological model should be determined from experimental and statistical approaches,
which will be conducted in the following study. Furthermore, in MINoSGI, daily photosynthetic rate was determined as a result of energy, water and carbon exchanges
between the atmosphere and forest as affected by variations in meteorological and hydrological conditions suchas light, temperature and soil water from hour to hour.
The phenological framework of the present study that anaccumulation of the assimilated carbon into the carbon
pool at the end of growing season is used to build leaves
only for the next year is not necessarily true of all species.
Hence, further implementation of the determination of
maximum LAI relevant to the seasonal cycle of carbohydrate mobilization in the tree body might be required in
the present phenological scheme, by considering speciesspecific features.
Numerical experiments confirmed the importance of thepresent canopy development model, i.e. determining theannual variation in LAI of a deciduous forest based on the
carbohydrate reserve and mobilization, and has enhanced
the possible utility of the forest dynamics model,
MINoSGI, coupled with the canopy development model.
Furthermore, with the fully validated model, we foundthat doubled CO2 concentration increased and increased
air temperature reduced the production of a deciduous
forest over the 150-year simulation period. Hence, thesensitivity test suggests that the deciduous forest ecosystem in the cool-temperate climate zone has enhanced efficiency of carbon absorption over a hundred years indoubled CO2 at increased temperatures by <5°C compared with the present climate conditions but deteriorated
one in doubled CO2 concentration at increased temperatures by ≥ 5°C. The former is due to the CO2 fertilization
effect that surpasses the increased temperature's negative
effect on production (i.e. respiration is increased by temperature), the latter vice versa.
As a next step, using the improved MINoSGI coupled withthe canopy development model with variable LAI asaffected by varying environments, we can predict variations in the annual carbon budget of the deciduous forestsin the cool-temperate and boreal climate zones, where
there is great apprehension about the effects of global
warming on forest ecosystems, with climate change scenarios. Therefore, we can expect not only to elucidate therole of deciduous forests of these regions in the global carbon cycle but also to apply our model to assessing of the

Page 7 of 18(page number not for citation purposes)

Carbon Balance and Management 2007, 2:6

effects of forest management and climate change on the
carbon stocks in managed and unmanaged (or natural)
forest ecosystems.

Methods
1. MINoSGIFor the purpose of describing the interactions between
meteorological physics and the growth of plants with theprogress of time, Watanabe et al. [28] developed a gap
model combined with a multilayered microclimate environment and forest dynamics, MINoSGI. MINoSGI consists of two submodels, i.e., a plant size dynamics modeland a micrometeorological model, and simulates dynamics such as the change of plant growth and mortality due
to competition between individual trees within a given
physical environment inside a forest. Such dynamics ofthe forest, in turn, cause changes in the physical environment. This interactive feedback process allows us to predict the impact of climatic change on forest dynamics on
a long time scale. Watanabe et al. [28] applied MINoSGI
to an evergreen coniferous Sugi stand (Cryptomeria
japonica D. Don, the Japanese cedar, Taxodiaceae) inmonoculture, and was successful in reproducing the sizestructure dynamics of the forest over six years in comparison with ecological inventory data obtained from 1983 to
1988 [42]. For the purpose of prognosticating physicalvariables such as soil temperature and moisture, snow
mass, snow surface temperature, snow albedo and runoff
in the soil ground process, MINoSGI includes the groundprocess scheme similar to MATSIRO (Minimal Advanced
Treatments of Surface Interaction and RunOff), which is aland surface model developed to be coupled with anatmospheric general circulation model [43]. In the
ground process scheme, five vertical multi-layered distributions of soil temperature and moisture are prognosticated in the soil for a whole forest stand. Snow surface
temperature and albedo on the ground are also predictedby this model. These are calculated from the input variables such as nonconvective rainfalls and snowfalls using
the fluxes in MINoSGI at hourly time steps. With prognosticated soil moisture, soil water stress was calculated inaccordance with a stress function due to soil water deficit
defined by [44]. Thus, the inclusion of soil water stress in
the ground process scheme into MINoSGI affects dailyleaf photosynthesis, namely, the height- or diametergrowth of trees under a specific soil moisture condition.
The treatment can be acceptable from the observational
result that the increment of diameter growth of trees isconstrained by soil moisture deficit [45].

The fundamental treatment to predict forest dynamics in
MINoSGI is described as follows. In MINoSGI, the prognostic equations predicting the dynamics of a distribution
density f(t,h) of individuals and the total plant mass

http://www.cbmjournal.com/content/2/1/6

f(t,h)w(t,h) of plant height class h per unit ground area attime t are given by∂f ( t , h )
∂t

=

1 ∂2
{ D ( t , h ) f ( t , h )} − ∂∂h { G ( t , h ) f ( t , h )} − M ( t , h ) f ( t , h ) ,
2 ∂h2

(1)
and
∂f ( t , h ) w ( t , h )∂t= Pn ( t , h ) f ( t , h )+1 ∂2{ D ( t , h ) f ( t , h ) w ( t , h )} − ∂∂h { G ( t , h ) f ( t , h ) w ( t , h )}2 ∂h2

− M ( t , h ) f ( t , h ) w ( t , h ) − f ( t , h ) ⎡⎣ L f ( t , h ) + Lw ( t , h ) + Lr ( t , h ) ⎤⎦ ,(2)
where w(t,h) is mean plant mass of each h, D(t,h), G(t,h)
and M(t,h) are variance and mean of growth rate inheight, and mortality rate of individuals, respectively.
Pn(t,h) is the net primary production (NPP) per plant ofeach plant height class and Lf(t,h), Lw(t,h) and Lr(t,h) are
the litter production rate for leaf, wood and fine roots,
respectively. These equations include the original diffusion model proposed by [29], which describes fluctuations caused by various factors such as spatial
heterogeneity of environmental variables, spatial variation of individuals in the neighborhood competitiveeffects [29,46]. However, in the present study, the diffusion term D(t,h) in eqns. (1) and (2) is excluded because
it is assumed that the experimental forest area is coveredby a spatially homogeneous land surface.
In MINoSGI, the vertical distribution of leaf area is estimated from the pipe-model theory proposed by [47],which relates leaf area of plant height h and stem diameterd of a plant. The vertical profile of stem diameter of plantheight h, d(t,z,h), is described as follows:⎡ ⎛ z ⎞η ⎤d ( t , z , h ) = d0 ( t , h ) ⎢ 1 − ⎜ ⎟ ⎥ ,⎢⎣ ⎝ h ⎠ ⎥⎦(3)

Where z is height from the ground, d0 (t,h) is stem diameter at the ground and η is a parameter. As η increases, theshape of stem changes from conic to spheroidal [48].
With the stem diameter distribution, the vertical distribution of leaf area is given by

FL ( t , z , h ) =( ) = 2θηd 2 ⎡⎢1 − ⎛ z ⎞η ⎤⎥ zη −1∂ θ d2
∂z

0

⎢⎣

⎜ ⎟
⎝ h ⎠ ⎥⎦ hη(4)with the allometric parameter θ and the vertical distribution of potential leaf area of a height class h, FL(t,z,h). In a

Page 8 of 18(page number not for citation purposes)

Carbon Balance and Management 2007, 2:6

http://www.cbmjournal.com/content/2/1/6

growing season

1
0
Ds

De

Ds

De

day of year

Stem diameter

Carbonreserve ratio

growing season

LAI(Leaf area index)

(iii)

(i)
(iv)

(ii)

(i)

(iii)

(ii)

day of year

LAI max(iv)

LAI max

tb

t0

tf

1st. yeartb

t0

2nd. yeartftb

day of year

3rd. year

Figure 6 diagram
Schematic
deciduoustreeof the seasonal courses of LAI, radial growth and the carbon reserve ratio with phenological events for a
Schematic diagram of the seasonal courses of LAI, radial growth and the carbon reserve ratio with phenological events for a
deciduous tree. t0, tb and tf are the onset of bud dormancy, the timing of bud burst and the completion of leaf fall, respectively.
(i) is the combined duration of autumn and winter dormancy. The duration of (ii) is the ontogenesis development stage toward
bud burst. In addition, the duration of (iii) and (iv) are the stages of LAI increment from bud burst to its maximum and LAI decrement from the maximum to the completion of leaf fall. The timings Ds and De mean the onset and end, respectively, of carbon
reservation in the storage pool of woody nonphotosynthetic organs. For simplicity, Ds was assumed to be identical to De (i.e.,
Ds = De = Dx) in this study.

dense forest, the actual leaf amount at the bottom of the
canopy layer is reduced by competition between individuals, i.e., self-pruning or the rising of crown. Watanabe etal. [28] calculated the actual leaf area FLA(t,z,h) by incorporating the temporal change of leaf expansion coefficientfor each layer at daily time steps based on a leaf carbon
budget. When a negative carbon budget per leaf of a specific layer is calculated, FLA(t,z,h) at the layer is estimatedby reducing the leaf amount of the corresponding dry
matter by the negative carbon amount. In the frameworkof this study, FL (t,z,h) is calculated by considering the difference between maximum tree height of a height classand the height of the mean crown base zcb, which is
decided as the height after the occurrence of the crown risedue to the performance of dead leaves. Using this FL (t,z,h)as the upper limited leaf area, the actual maximum leaf
area FLA(t,z,h) is determined through the comparison withleaf area, which is produced by phenological processes inwhich the leaf area is determined by using carbon assimilation by tree leaves stored during the latter part of thegrowing season until defoliation to build a leaf amountfor the next year.
The model parameters, which cannot be determinedbecause of the lack of observational approaches and/or
are undocumented, are adjusted into formulated relationships using inventory data. The vertical profile of crowndiameter dcrown is calculated as
1

z⎞ 2
⎛
(5)
dcrown ( t , z , h ) = δ d0 ⎜ 1 − ⎟h⎠⎝with the adjusted parameter δ, stem diameter at the
ground d0 and distance from the ground level z. The next
adjusted parameter fd is used for the relationship between

Page 9 of 18(page number not for citation purposes)

Carbon Balance and Management 2007, 2:6

http://www.cbmjournal.com/content/2/1/6

leaf respiration rate Rd and maximum catalytic activity of
Rubisco Vm as follows [28]Rd = fdVm.

(6)

In addition, we considered the allocation relationshipbetween plant height h and stem diameter d0 in eachheight class. The difference in allocation between individual tree height classes might produce different competition modes, resulting in different forest dynamics.
Watanabe et al. [28] introduced the following relationshipon the basis of [49]:

⎛
⎞
Δh
1
= μ0 + μ1 ⎜
−1⎟,
⎜
⎟
Δd0
⎝ Im ( t , h )⎠(7)where μ0 and μ1 are allocation parameters, and Im(t,h) is
the simulated value of the absorbed relative photosynthetically active radiation (PAR) intensity averaged overthe entire leaf area of a single plant of height h at time t.
From eqn. (7), it is apparent that the allocation relationship is dependent on the light environment within a forest stand. Hence, a larger μ1 might enhance the height
growth of smaller plants in a dark light environment
inside the forest, consequently reducing differences ingrowth rate between trees of different height classes. Thus,
these parameters can be used as a measure of competition
for light between trees [49]. The general outline andmodel components of MINoSGI are described in detail by
Watanabe et al. [28].
2. Soil carbon dynamics
The net carbon uptake in a forest ecosystem, net ecosystem production (NEP), is given as

NEP = (GPP - Ra) - Rh = NPP - Rh,(8)

where Ra and Rh are autotrophic and heterotrophic respiration rates, respectively. Ra is a carbon loss by a plant, thesum of growth respiration and maintenance respiration ofnonphotosynthetic (sapwood and fine roots) organs [33].
In contrast, Rh is the loss of carbon accompanying with theactivity of soil microorganisms that decompose soilorganic matter (SOM). In the original MINoSGI, Rh wasgiven as a function of soil temperature without considering carbon cycling in the soil [28]. Here, we incorporated
a below-ground soil carbon dynamics model into
MINoSGI and Rh is now modelled through the decomposition processes of soil organic carbon related with theabove-ground leaf fall event. The framework of the soil
carbon dynamics is based on [2]. Litter and two SOM poolcompartments, each of which has a different decomposition rate, were considered. Decomposition of each pool iscalculated as a function of soil temperature and moisture
at hourly time steps. A fraction of decomposed litter isreleased to the atmosphere and the remaining enters the
SOM pools. The allocation rate to each SOM pool is given
as a constant. Then, Rh is obtained as the summation of
respiration rates from each compartment. All parameters
used in the soil carbon dynamics model are the same as in[2].
3. Description of the deciduous canopy development
model
The scheme of processes incorporated into the present
deciduous canopy development model is shown in Fig. 6.
The canopy development model is composed of a phenological model and a model for maximum leaf area index
based on the carbon reserves in woody tissues.

The phenological model describes the process of chilling
requirement during dormancy that must be fulfilled
before the ontogenetic development toward bud burstcommences, and determines the timing of bud burst
through the stage of ontogenetic development after the
release of dormancy [50]. In Fig. 6, tb and tf indicate thedates of bud burst and completion of leaf-fall, respectively. Duration (i) is the development stages of bud dor-Table 1: The parameter list for the Larix kaempferi stand.

SymbolDefinition

Value

Unit

Reference or adjusted in the present study

Vm
α
θ
SLA
η
δ
fd
μ0
μ1
Dx
Dcrit
Bcrit

Maximum catalytic activity of Rubisco
Allometric parameter
Allometric parameter
Specific leaf area
Allometric parameter
Factor of crown diameter
Ratio of leaf respiration to Vm
Parameter for light competition
Parameter for light competition
Parameter for the start of carbon reserves

35
0.022
0.0256
16.7
1.9
4.0
0.017
0.3
0.5
230
50
80

μmol m-2 s-1
kg cm-2 m-1
m2 cm-2
m2 kg-1
m cm-1-

Kurachi (1986)
Kurachi (1986)
Kurachi (1986)
Kurachi (1986)
adjusted in eqn.(3), (4)adjusted in eqn.(5)adjusted in eqn.(6)adjusted in eqn.(7)adjusted in eqn.(7)Parameter for threshold value to determine the date of bud burst

This issue (fixed)
This issue (fixed)

Page 10 of 18(page number not for citation purposes)

Carbon Balance and Management 2007, 2:6

http://www.cbmjournal.com/content/2/1/6

mancy period (autumn-winter dormancy) starting at t0(thus buds growing in the next year are assumed to beformed before t0), where we define 'dormancy' as a statein which buds do not grow even under favorable ambientconditions. Normally, the dormancy processes between
autumn and winter are different. Thus, each has a different response function of temperature [15,51,52]. However, in recent studies, the same developmental response
function of temperature has been used consistently
through the autumn to winter periods [15]. (ii) is a stage
during the ontogenetic development of buds toward bud
burst after duration (i) has been completed. (iii) and (iv)
are the increment and decrement processes of the leaf
amount in a year, respectively. Ds represents the date on
which trees start to store part of surplus assimilates in thetree body, and De the date after which all the surplus
assimilates are stored without being used for tree growth.
In the present study, we assumed to be Ds = De≡Dx andtreated Dx as a fixed date. The outline of the above-mentioned processes is specified in the following. Details on
the treatment of Ds and De are described later.
To determine the timing of bud burst, approaches proposed by [51,52] were used: these include the determination of the stage of dormancy to describe the chillingrequiring stage indispensable for the release fromautumn-winter dormancy, and the stage of ontogenesis
essential for the induction of bud burst after the autumnwinter dormancy period.
The time change of the stage of dormancy SD(t) can be
written in the form aggregated with time t using the rate ofbud development fD(t) as follows:t

SD ( t ) = ∫ fD ( t ) dt ,t0

(9)where t0 is the beginning of the autumn-winter dormancy.
Häkkinen, Linkosalo and Hari [14] applied the dependence of fD(t) on temperature T(t) based on experiments by[52] to seedlings and seeds of birch (Betula pubescens), andfD(t) = fD(T(t)) was described as⎧0⎪
T ( t ) − TmfD ( T ( t ) ) = ⎪⎪⎨ To − Tm⎪⎪ T ( t ) − TM⎪⎩ To − TMT ( t ) ≤ Tm or T ( t ) ≥ TM
T ( t ) > Tm
T ( t ) < TM(10)
with the optimal chilling temperature, To (To = 6.9), while
Tm and TM are constants (Tm = -3.4, TM = 10.4) [17]. On theother hand, the stage of ontogenesis So(t) is described by

the rate of ontogenetic development, fo(t) based on [51]([14]) as follows:t

So ( t ) = ∫ fo ( t ) dt ,

(11)

tf

where tf is the beginning of the stage of ontogenesis. In
addition, fo(t) is given as
fo(t) = go(T(t))·wo(L(t)),(12)

where go(T(t)) denotes the dependence of the rate of
ontogenesis on temperature as follows [51,17]

⎧=
0
⎪
28.4
go ( T ( t ) ) ⎨=⎪−0.185( T ( t )−18.4 )1+ e⎩

T(t ) > 0T(t ) > 0,
(13)

and wo(L(t)) means the effect of the light conditions (suchas light intensity, night length, and spectral compositionof light) on ontogenesis. L(t) denotes the signal from light
conditions [14]. Changes in light conditions during theannual cycle are the most reliable source of informationreflecting the timing of seasons, especially in the boreal
zone because these seasonal changes are much larger thanthose in lower latitude areas and the photoperiod plays arole in controlling the timing of phenological events ofvarious boreal tree species [14,15]. The effect of thesechanges was considered in the present case of cool-temperate climate area, i.e. we assumed that wo(L(t)) is givenby⎪⎧ = 0 L ( t ) < Swt(14),
wo ( L ( t ) ) ⎨⎩⎪ = 1 L ( t ) > Swt
where Swt is daily mean solar radiation. Häkkinen, Linkosalo and Hari [14] used calendar date as an operationalvariable for light conditions. Instead of calendar date, weused Swt and assumed to be wo(L(t)) = 1 when Swt was
greater than 300 Wm-2 which can be obtained for theperiod from the middle of February to the end of Octoberin the present study site under sunny weather conditions.
Among the above functions, the unknown variables necessary to be established are the dates of the beginning andend of the autumn-winter dormancy, and a decisionneeds to be made as to when ontogenetic development ofbuds toward bud burst is completed. Normally, the onset
of autumn-winter dormancy is likely to vary with species
and annual environmental variables, and should therefore be estimated using an experimental approach [51] orby statistical interpretation using long-term phenological
data [14]. In this article, it is assumed that the beginning
Page 11 of 18(page number not for citation purposes)

Carbon Balance and Management 2007, 2:6of the autumn-winter dormancy in the process (i) in Fig.
6 occurs on the same fixed date each year as in [52]. Whilethe end of the autumn-winter dormancy is describedusing the accumulated sum of the value of SD(t) showing
the state of dormancy obtained from the response function of air temperature, the stage is completed when its
accumulated sum surpasses an assigned threshold value(Dcrit). The process in (ii) is also determined by the sum ofvalues showing the developmental process of the bud
toward bud burst, calculated by the logistic response function of air temperature just after the dormancy period (i).
Bud burst occurs when S0(t) in Eqns. (11) exceeds a given
threshold value (Bcrit). These threshold values assigned in
Eqns. (9) and (11) to determine the dates of the end of theautumn-winter dormancy and bud burst are establishedbased on the temperature-dependence of leaf opening
obtained from inventory data in Japanese larch site (L.
kaempferi spp.) by [53] (see Table 1).
In addition, the incremental process (iii) of leaf amountafter bud burst until LAI is at a maximum change with thethermal condition in the current year. We use the growing
degree-day as an index of the thermal condition which isdefined as the summation of the canopy-effective temperature (Tce), where Tce is given by the subtraction of the basetemperature (herein, the base temperature is set to 5°C)
from the canopy-averaged temperature. The canopy-averaged temperature is estimated by dividing the summationof a product of leaf amount and temperature at each vertical layer obtained from multilayered meteorologicalmodel by total leaf amount over a canopy. When the accumulated temperature is greater than an assigned thresholdvalue which is set to 250°C [54,55], LAI reaches a maximum. On the other hand, leaf-fall in autumn throughduration (iv) is described using a 10-day running mean
temperature calculated from the daily mean temperature,and the start of leaf-fall is determined when the 10-day
value is less than an assumed threshold value, and allleaves fall within 30 days.
4. Contribution of carbon reserve to the maximum LAI- a
maximum LAI model
The leaf amount of a tree has a close allometric relationship with the stem diameter at the bottom of the crownheight (the pipe model; [47]) as denoted in eqn.(4). However, it is likely that the leaf amount of a deciduous tree is
not only related to the stem diameter but also closely tothe carbohydrate reserve in the trees (e.g. [56,57]). Determination of the maximum leaf amount (i.e., maximum
LAI) is greatly affected by both the environmental variables of the current year and the carbon storage of the previous year which in turn reflects the environmental
conditions of the previous year. Therefore, the decisionsof Ds and De in Fig. 6 (when part of surplus production
starts being stored in the tree body, and when all surplus

http://www.cbmjournal.com/content/2/1/6

production starts being stored) allow us to make explicit
the effect of environmental differences on the annual
amount of carbon storage.
The instant that bud burst occurs, deciduous trees startphotosynthesizing with their unfolded leaves in the current year. The increase of leaf amount in spring, however,is likely to be so rapid that photosynthesis by leaves is notenough to support the rapid growth rate [45]. Hozumi
and Kurachi [57] investigated the relationship between
the relative growth rate of opened leaves and the seasonality of leaf biomass production using a component
model approach, and suggested that the total growth rateof leaf biomass in spring was largely supported by the carbon supply from non-photosynthetic organs where carbon in the previous year is stored. Moreover, Hozumi and
Kurachi [57] found that leaf amount was increased due toboth the contributions of carbohydrates reserved in theprevious year and net assimilates of unfolded leaves in thecurrent year through the seasonal period until the leaf
amount is at maximum.
To combine the effect of carbon stores and the pipe model
theory as concomitant factors determining the leaf
amount of a tree, the following was assumed. If the carbohydrates are abundant (i.e., if the dry weight of carbohydrates reserved in the tree body is greater than the leaf dry
weight of the tree predicted from the stem diameter
through the allometric relationship of the pipe model),the maximum possible leaf amount of a tree is equivalent
to the predicted value. But if the carbohydrates are notenough, realized leaf amount of a tree is smaller than the
one predicted by the pipe model. It should be noted herethat the dates of Ds and De in Fig. 6 is important to explicitly include the effects of environmental variations on theamount of carbon storage, because these dates determine
the length of duration during which surplus production is
accumulated as carbon stores. Hence realized LAI of theforest may also be affected by the carbon storage of thetrees in the previous year, which in turn shows that theenvironmental conditions of the previous year affects thestand LAI in the current year. Actually, LAI of a forest canvary from year to year with climatic conditions and environmental events.
As mentioned above, it is incidentally known that part ofthe carbon resources is also used for spring woody production [45] and reproductive organs. Moreover, thestored carbohydrate (starch) might be used for active
radial (or cambial) division (new cell-wall formation)
before the full development of leaves after bud burst [58].
However, we ignored such processes in this study to simplify the relationship between the leaf amount and thecarbon resources. The assumption that leaf growth comes
only from the carbohydrate reserve while plant growthPage 12 of 18(page number not for citation purposes)

Carbon Balance and Management 2007, 2:6

comes only from newly assimilated carbon might be difficult to apply to all tree species. For instance, it might betrue of diffuse porous trees like birch, whilst not necessarily true of ring porous trees like oak where spring wood
growth must start before bud burst to support leaf growth[45]. The observational result on the translocation fromcarbon reserves in spring of L. kaempferi species by [57]might support the treatment on the carbon translocationas part of ring porous trees, although L. kaempferi does notbelong to them. Nevertheless, in this study, we employed
this assumption to describe explicitly the possible relationship between maximum LAI and carbon reserves.
Aber and Federer [19] developed a forest carbon and waterbalance model (PnET) to describe monthly time-steps ofwater and carbon balances in a forest, and Aber et al. [18]further improved the original model (PnET-II), in which
bud and wood carbon storage pools accumulated on amonthly basis during several years before the current year
were used to simulate the foliar and wood production for
the next year, respectively. Maximum leaf amount (i.e.
LAI) for the next year was then decided, and PnET-II withthe routine of carbon storage was applied to hardwood,
spruce-fir, mixed hardwood/spruce-fir and mixed hardwood/pine forest sites [18]. In the present study, we used
L. kaempferi which has short and long branches. In general, leaf primordium of L. kaempferi species during the
middle of growing season in a given year is likely to perform on a short branch for leaf production in the following year, and the formation of leaf primordium of Larix
tree species such as L. kaempferi and Larix gmelini appears
to be largely affected by varying meteorological conditions before the onset of the formation (personal communication with Dr. Uchiyama from her unpublished data).
After the leaf expansion in the next year, most of the carbon assimilation by opened leaves on short branches isprovided for the growth of long branches in the corresponding year. Additionally, leaf primordium is formedon long branches, and new leaves from the long branches
occur after all leaves from leaf primordium on shortbranches are opened. It implies that most of carbon
assimilation by the leaves in the current year might beused to construct additional buds and wood productionin the year along with the rest of carbon reserves stored in
the previous year. However, it is difficult to separate how
carbon amount from reserves in a previous year andgained by opened leaves in the current year to produce
wood and foliar growth has been used. In our study, it was
assumed that the production of leaves in the current year
depended on the carbohydrates from the previous year
which contribute to the entire leaf amount (i.e. maximum
LAI) over a stand, including leaf amount from the secondary wood production (long branches) produced by theopened leaves on short branches. Further, we assumed
that carbon assimilation by leaves produced in spring is

http://www.cbmjournal.com/content/2/1/6

used for plant buds and wood growth in the current year
after bud burst. With these simple treatments andassumptions to link the woody carbon storage with the
maximum LAI, we modeled phenological processes ofdeciduous L. kaempferi.
5. Modeling carbon storage in association with seasonalcarbohydrate mobilization
After the leaf amount reached a maximum, part of theassimilated carbon may continue to be utilized in further
growth. The rest of the carbon may possibly be transformed into carbohydrates such as sugar (including glucose, fructose and sucrose), and starts to be stored as the
carbon pool in the woody tissues and bark of branches,trunks and roots during the periods of low photosynthesis
to fuel the maintenance of respiration, especially during
autumn [59]. These surplus photosynthetic products are
stored until the end of the growing season [57]. Seasonalreserve mobilization of carbohydrates has been describedin several types of trees and for different stand ages
[45,56] and the seasonality has been quantitatively investigated (e.g., Betula ermanii by [60]). Hence we assumed
that reserved carbohydrates are used for respiration oftrees in leafless seasons and construction of leaves in thenext spring.

The total amount of carbon stored in the tree body isaffected by the timing of growth cessation of trees, as
stated before. Stem-diameter growth (radial growth) often
continues for a long time after shoot elongation ceases inlate summer to autumn [61]. In addition, the seasonalpattern of shoot elongation and diameter growth in B.
ermanii was observed during the growing season over two
years (2003–2004) in Hokkaido, north of Japan, and it
was found that cessation of shoot elongation was earlier
than that of radial growth (unpublished data by M. Toda,
K. Kato and K. Laska). Moreover, the timing of diameter
growth cessation was dependent on tree size (tree height),
but it mostly occurred around the beginning of September
in both years. The tendency was consistent with the reportby [53], which showed that diameter growth in L. kaempferi stopped around the end of August or the beginning of
September. Thus, in the present model, it was assumed
that all surplus assimilated carbon after cessation of diameter growth is stored in the body of the tree. As for thedates of Ds and De, however, we assumed in the present
paper that all the surplus production starts being stored
on a certain day of year (i.e., Ds = De≡Dx in Fig. 6), andtreated Dx as a fixed date. For estimating Dx, we simply
referred to the end of diameter growth based on the observational fact that the diameter growth of L. kaempferi stops
around the end of August or the beginning of September[53]. We treated Dx as an adjusted parameter in this study
and determined it as follows. For a given range of Dx, thesize distribution of plant height class was simulated

Page 13 of 18(page number not for citation purposes)

Carbon Balance and Management 2007, 2:6

http://www.cbmjournal.com/content/2/1/6

together with the adjusted parameters in eqn.(5), (6), (7);Dx was then optimized by finding Dx that gives the best
agreement between the observed and simulated size distributions.

'84 observed'84 simulated

PAR/PAR0

1

0.1

0

2

4

6

8

Cumulative leaf area (m2m-2)
Figurecumulative
feri
stand7 inand
leaf
1984
area perunit
groundarea
in relative
the 0Larix) against
kaempObservedofphotosyntheticallysimulated
activeverticalradiation
profiles(PAR/PARofintensity
Observed and simulated vertical profiles of relative intensityof photosynthetically active radiation (PAR/PAR0) against
cumulative leaf area per unit ground area in the Larix kaempferi stand in 1984. PAR in each layer is normalized by PAR at
the canopy top (PAR0). Six observed data were obtained
within the stand in each year.

Tree density (m-2)

0.4

observed
simulated

6. Ecological inventory and meteorological data sets
Ecological inventory data from 1982 to 1989 used for themodel validation were obtained from a plantation of Japanese larch (Larix kaempferi (Lamb.) Carrière), a deciduous tree species, in the Nagoya University Experimental
Forest at Inabu, Aichi Prefecture in central Japan. Total
precipitation during the growing season from April to
October was 1494 mm and the mean air temperature was
7.6°C over the three years from 1982 to 1984 [62]. The
study plot of 200 m2 (20 m × 10 m, 35°11'N, 137°33'E,
1040 m a.s.l.) was established in 1982 [63]. In 1982,
stand age was 18 years. Tree density was 2580, 2529 and2327 trees ha-1 in 1982, 1983 and 1984, respectively. The
mean height of the trees, the mean stem diameter at breast
height, and the mean stem diameter at the lowest branchheight increased from 10.3 m to 11.7 m, 10.7 cm to 11.8
cm, 5.2 cm to 6.0 cm during the three years, respectively.
More details are described in [62].

Meteorological data are required as a set of forcing data
for MINoSGI. Unfortunately, simultaneous observation
procedures were not established with the above-mentioned ecological inventory data. Hence, use was made ofa specified data source from the AMeDAS (AutomatedMeteorological Data Acquisition System) at the station at
Inabu (35°12.6'N, 137°30.6'E, 505 m a.s.l.) 5.7 km awayfrom the site. The same data correction for altitude as usedby [28] was applied to obtain air temperature and pressure with altitude assuming a moist adiabatic lapse rate of
0.005 Km-1, and humidity calculated from the correctedvariables. Because the present study site is very close to thesite used by [28], the other meteorological variables usedin the present study were the same as those described by
[28].

0.3

0.2

0.1
1982 1983 1984 1985 1986 1987 1988 1989 1990

Time (year)
FigureObservedLarixkaempferi
8 and simulated
stand fromannual1982changes
to 1989in tree density of the
Observed and simulated annual changes in tree density of theLarix kaempferi stand from 1982 to 1989.

7. Simulation setup
The inventory data include the number of trees, the treeheight (h) and stem diameter at breast height (DBH) andat the lowest branch height (DB) and the total weight (w)for an individual tree. Using time steps of one year for thesize dynamics model and one hour for the microclimate
model, simulation over seven years from 1982 to 1989
was carried out. Thus, structural changes of the forest were
derived at yearly intervals. However, leaf photosynthesis,
calculated on the basis of a combination of the photosynthesis at the biochemical level [64], the stomatal conductance of the semi-empirical model proposed by [65], therespiration at a leaf scale in eqn. (6) together with thewoody and root respiration in the microclimate model,
can be calculated hourly for an area-averaged and multi-

Page 14 of 18(page number not for citation purposes)

Carbon Balance and Management 2007, 2:6

Frequency (m-2)

0.08(a) 1983

1982 observed
1983 observed
Simulated

http://www.cbmjournal.com/content/2/1/6

0.08(b) 1985

1985 observed
Simulated

0.08

0.06

0.06

0.06

0.04

0.04

0.04

0.02

0.02

0.02

0

0

0

0

5

10

15

20

25

0

5 10 15 20 25
Tree height, h (m)(c) 1989

0

5101989 observed
Simulated

15

20

25

Figure 9 and simulated frequency distributions of Larix kaempferi trees in each height class in (a) 1983, (b) 1985 and (c) 1989
Observed
Observed and simulated frequency distributions of Larix kaempferi trees in each height class in (a) 1983, (b) 1985 and (c) 1989.
All live trees were classified in tree height at 1-m intervals. The vertical axis denotes the number of trees in each tree height
per study area (200 m2). In Figure 9 (a), the initial size distribution in 1982 is also presented.

layered space. Therefore, the framework enables us to estimate the seasonal production at daily intervals, whichleads to the determination of phenological events in thespring of the subsequent year as specified in Methods. The
inventory data from 1982 were used as the initial input
data, and seven-year simulated outputs were comparedwith the observational data for each year. The simulation
model was given boundary conditions of tree height classwith zero flux at the lowest boundary and no outflow flux
at the top boundary, during the simulation period. Theseliving trees were grouped at 1-m intervals into heightclasses from 0 m to 30 m for each year. In order to produce
a proper calculation, the closure of annual carbon balancehas been examined for each one-year time interval.
The species-specific model parameters in MINoSGI
obtained from the inventory data in the study plot arelisted in Table 1. In addition, the adjusted parameters arealso summarized. In the MINoSGI simulation, we used
the same model parameters for respiration of stem andfine roots of L. kaempferi as those of C. japonica used in theprevious research [28] because of the lack of inventory
data for these parameters. Hukami [66] and Hagihara andHozumi [67] investigated the difference in the annual respiration rates of stem and fine roots between L. kaempferi
and C. obtusa (which is one of typical coniferous species
and has eco-physiological characteristics similar to C.
japonica) at the stand level. Hukami showed that stem andfine root respiration rates of 16- to 24-year-old L. kaempferi were 4.2–8.4 and 4.3–5.4 ton ha-1 year-1, respectively.
Hagihara and Hozumi [67] showed that those of 18-yearold C. obtusa were 9.9 and 4.6 ton ha-1 year-1. Therefore,

we assumed that the difference in stem respiration rate
between these species does not largely affect the production variables over seven years during the present simulation period because stand age of L. kaempferi during theseven years is 18 to 25 years. All of the model parametersof except for those given in Table 1 were therefore thesame as in [28]. Furthermore, we determined initial carbon stocks in litter and two SOM pools in soil carbon
dynamics model based on the soil respiration measurement in C. japonica stand in the Inabu experimental field[67]. Initial litter stock was 3.1 tC ha-1. On the other hand,values of the reminder initial SOM C stocks were given as
6.0 tC ha-1 in intermediate and 7.6 tC ha-1 in slow SOM
pools, respectively, from estimation of mineral soil androot components by [67].
8. Determination of model parameters and validation
Some unknown parameter values were adjusted by simulation because these were difficult to obtain from observational approaches or due to a lack of observations. Of theunknown parameters described in Table 1, the four ones
treated in [28] were investigated; δ in eqn. (5), fd in eqn.
(6), and μ0 and μ1 in eqn. (7) (see Methods section). The
value of η in eqns. (3) and (4) was also treated as one of
the adjusted parameters. In addition, as described in
Methods section, the additional parameter that characterizes deciduous trees in conjunction with phenologicalevents was simultaneously adjusted with the above-mentioned variables.

For instance, the parameter fd affects the time change ofstand-level variables such as tree density, averaged tree

Page 15 of 18(page number not for citation purposes)

Carbon Balance and Management 2007, 2:6

height, LAI, net primary production (NPP), biomass andsize distribution of tree height. Larger fd brings about anincrease in leaf respiration and therefore may increase themortality of smaller-sized trees, which generally gain
lower carbon assimilation due to the low light intensity
they receive in the forest. On the other hand, the timechange of size distribution and LAI did not always depend
on fd alone. The simulated size distribution is largely
dependent on the height-diameter growth relationshipparameters (μ0 and μ1), which change the light environment within the forest. In addition, the timings of the startof carbon reservation (Dx) may have an effect on LAI andsize distribution. Thus, an interactive adjustment between
these parameters should be made for their final estimation by comparing the simulated and observed results ofvertical profile of light intensity within a canopy (Fig. 7),
tree density (Fig. 8), LAI (Fig. 2(a)) and size distribution(Fig. 9).
Initially, we assumed that the crown of L. kaempferi is parabolic-shaped. The parameters δ, η were determinedwhen the simulated vertical profile of light intensity withdepth into the canopy is consistent with the hourly-averaged measurement on photosynthetically active radiation(PAR) observed on specific days in 1984 by [34]. In thisoperation, we gave a constant allocation parameter set of
μ0 = 0.5, μ1 = 0.5, fd = 0.017 and Dx = 200. μ0 and μ1
may be involved in the change in light intensity within acanopy, but their effect on light intensity was very small
compared to those of δ and η. When the combination ofδ = 4.0 m cm-1 and η = 1.9 was chosen, the simulated andobserved values were highly correlated with R2 = 0.83.
The magnitude of parameter δ represents the ratio of theground-level crown diameter to the ground-level stem
diameter for an isolated plant [28], and we obtained aratio of 40. The value of η = 1.9 represents the canopy
shape with widely vertical-attached leaves from the bottom to the top of an isolated tree [48], and the crownshape agreed with the vertical shape of a L. kaempferi tree.
Furthermore, fd was adjusted so that the simulations
would reproduce the annual live and dead tree numbers
well (Fig. 8), both of which are the consequence of competition between trees with the progress of time. The valueof fd was estimated by parameter adjustment as fd = 0.017,
which is similar to that of C. japonica (fd = 0.022) [28]. In
this operation, we also gave a constant allocation parameter set of μ0 = 0.5, μ1 = 0.5 and Dx = 200. Simultaneously,with the parameter set, good agreement between theobserved and simulated maximum LAI in 1982 and 1984
was obtained (Fig. 2(a)) (at a particular time during thelate growth period each year; LAI observations of individual trees were carried out by [34]). Kurachi [34] estimatedthe value of LAI by using observed values of leaf biomass
and specific leaf area (SLA) in each vertical layer at 1-m

http://www.cbmjournal.com/content/2/1/6

intervals. Unfortunately, Kurachi [34] did not measure
SLA in upper layers near the top of the canopy in 1983,
and used SLA values measured in the lower layers to complement SLA in the upper layers. However, we found thatthe complemented SLA values in the upper layers in 1983
were smaller than the observed upper SLA values in 1982
and 1984. Therefore, we believe that LAI in 1983 obtainedby [34] was underestimated because the underestimated
SLA values were used (Fig. 2(a)). We assumed that theadditional parameter Dx should be indispensable to get a
final determination of fd. However, we found that the
impact of the difference in fd was larger compared withthat in Dx. Then, Dx was adjusted to DOY 200 (day of year)
for every year, and the observed and simulated dates ofbud burst over 3 years from 1982 to 1984 were 12 April
and 12 April, 26 April and 6 May, 15 April and 14 April,respectively.
Finally, we adjusted tuning parameters μ0 and μ1 so thatthe simulated size distribution of tree height would beconsistent with the observed results. As shown in eqn. (7),greater μ0 brings about a greater increase in tree growth in
the vertical direction compared with that towards the lateral diameter direction, whilst a greater value of μ1 a
greater height increment in the smaller-sized tree growth
in lower light intensity within the forest. The values of μ0
and μ1 obtained in this study were 0.3 m cm-1 and 0.5 m
cm-1. μ1 was extremely large in comparison with that of C.
japonica obtained by [28]. Although the estimation of theparameters was simply made so as to achieve agreement
between the observed and simulated results, differences inthe same parameter between tree species may be related tothe plant functional type, e.g., shade-intolerant species
such as deciduous L. kaempferi may have enhanced height
growth when the tree is small under a canopy, whereas
allocation to height growth in C. japonica, which is shadetolerant, may be less pronounced. The difference in theadjusted parameter fd may also affect the size distribution.
However, we found that the difference in fd does not
greatly change the size distribution.

Competing interests
The author(s) declare that they have no competing interests.

Authors' contributions
MT carried out numerical experiments and drafted themanuscript. MT, MY and TH conceived of this study, participated in its design and coordination; AS and TW
helped to draft the manuscript. All authors read andapproved the final manuscript.

Acknowledgements
We thank N. Kurachi for providing us with the inventory data used in thisarticle. We would also like to express our appreciation to K. Ono(Hokkaido Univ.) and T. Koike (Hokkaido Univ.) for helpful comments on

Page 16 of 18(page number not for citation purposes)

Carbon Balance and Management 2007, 2:6

this paper. Partial financial support was given by the 21st Century Centerof Excellence Program (E-01) funded by the Ministry of Education, Culture,
Sports, Science and Technology, Japan. This research was supported andfinanced, in part, by the Grant-in-Aid for Scientific Research to one of theauthors (MT). MT also received financial support from the Sumitomo Foundation, the Inamori Foundation, and the young scientist program of theInstitute of Low Temperature Science, Hokkaido University. We would liketo thank them all for supporting the research.

References
1.

2.

3.
4.
5.
6.

7.
8.
9.

10.
11.
12.
13.
14.
15.
16.
17.
18.

19.

20.

Pielke RA, Avissar R Sr, Raupach M, Dolman AJ, Zeng X, Denning AS:
Interactions between the atmosphere and terrestrial ecosystems: influence on weather and climate. Global Change Biology
1998, 4:461-475.
Sitch S, Smith B, Prentice IC, Arneth A, Bondeau A, Cramer W, Kaplan JO, Levis S, Lucht W, Sykes MT, Thonicke K, Venevsky S: Evaluation of ecosystem dynamics, plant geography andterrestrial carbon cycling in the LPJ dynamic global vegetation model. Global Change Biol 2003, 9:161-185.
Tans PP, Fung IY, Takahashi T: Observational constraints on the
global atmospheric CO2 budget. Science 1990, 247:1431-1437.
Baldocchi DD: The carbon cycle under stress. Nature 2005,
437:483-484.
Prentice IC, Heimann M, Sitch S: The carbon balance of the terrestrial biosphere: Ecosystem models and atmosphericobservations. Ecol Appl 2000, 10:1553-1573.
Zaehle S, Sitch S, Smith B, Hatterman F: Effects of parameter
uncertainties on the modeling of terrestrial biosphere
dynamics. Global Biogeochem Cycles 2005, 19:. GB3020, doi10.1029/
2004GB002395
Goulden ML, Munger JW, Fan SM, Daube BC, Wofsy SC: Exchangeof carbon dioxide by a deciduous forest: response to interannual climate variability. Science 1996, 271:1576-1578.
Arora VK, Boer GJ: A parameterization of leaf phenology forthe terrestrial ecosystem component of climate models. Global Change Biol 2005, 11:39-59.
Kramer K, Mohren MJ: Sensitivity of FORGRO to climaticchange scenarios: A case study of Betula pubescens, Fagus sylvatica and Quercus robur in the Netherlands. Climatic Change
1996, 34:231-237.
Rathcke B, Lacey EP: Phenological patterns of terrestrial plants.
Ann Rev Ecol Systematics 1985, 16:179-214.
Leinonen I, Kramer K: Applications of phenological models to
predict the futre carbon sequestration potential of boreal
forests. Climatic Change 2002, 5:99-113.
Matsumoto K, Ohta T, Irasawa M, Nakamura T: Climate change
and extension of the Gingko biloba L. growing season in
Japan. Global Change Biol 2003, 9:1634-1642.
Richardson AD, Bailey AS, Denny EG, Martin CW, O'keefe J: Phenology of a northern hardwood forest canopy. Global Change
Biol 2006, 12:1174-1188.
Häkkinen R, Linkosalo T, Hari P: Effects of dormancy and environmental factors on timing of bud burst in Betula pendula.
Tree physiol 1988, 18:707-712.
Leinonen I, Kramer K: Applications of phenological models to
predict the future carbon sequestration potential of boreal
forests. Clim change 2002, 5:99-113.
Kikuzawa K: A cost-benefit analysis of leaf habit and leaf longevity of trees and their geographical pattern. The American
Naturalist 1991, 138:1250-1263.
Chuine I, Kramer K, Hänninen H: Plant development models. In
Phenology: An Integrative Environmental Science Edited by: Schwartz MD.
Kluwer Academic Publishers; 2003.
Aber JD, Ollinger SV, Federer CA, Reich PB, Goulden ML, Kicklighter
DW, Melillo JM, Lathrop RG Jr: Predicting the effects of climatechange on water yield and forest production in the northeastern United States. Clim Res 1995, 5:207-222.
Aber JD, Federer CA: A generalized, lumped-parameter modelof photosynthesis, evaporation and net primary productionin temperate and boreal forest ecosystems. Oecologia 1992,
92:463-474.
Sellers PJ, Berry JA, Collatz GJ, Field CB, Hall FG: Canopy reflectance, photosynthesis, and transpiration. III. A reanalysis

http://www.cbmjournal.com/content/2/1/6

21.

22.

23.

24.
25.

26.
27.

28.

29.
30.
31.
32.
33.
34.
35.
36.

37.

38.

39.
40.
41.

42.

using improved leaf models and new canopy integration
scheme. Remote Sens Env 1992, 42:187-216.
Sellers PJ, Los SO, Tucker CJ, Justice CO, Dazlich DA, Collatz GJ,
Randall DA: A revised land surface parameteization (SiB2) for
atmospheric GCMs. Part II: The generation of global fields ofterrestrial biophysical parameters from satellite data. J Clim
1996, 9:706-737.
Kim W, Kanae S, Agata Y, Oki T: Simulation of potential impactsof land use/cover changes on surface water fluxes in the
Chaophraya river basin, Thailand. J Geophys Res 2005, 110:.
D08110, doi:10.1029/2004JD004825
Forely JA, Prentice IC, Pamankutty N, Levis S, Pollard D, Sitch S, Haxeltine A: An integrated biosphere model of land surface processes, terrestrial carbon balance, and vegetation dynamics.
Global Biogeochem Cycles 1996, 10:603-628.
Sato H, Itoh A, Kohyama T: SEIB-DGVM: A new dynamics global
vegetation model using a spatially explicit individual-based
approach. Ecol Model 2007, 200:279-307.
Ito A, Oikawa T: A simulation model of the carbon cycle inland ecosystems (Sim-CYCLE): A description based on drymatter production theory and plot-scale validation. Ecol
Model 2002, 151:143-176.
Kuroiwa S: Dry matter production of plants, Ecology and Evolution. In A Series of Modern Biology Iwanami Shoten; 1996. (In Japanese)
Barbaroux C, Breda N, Dufrene E: Distribution of above-groundand below-ground carbohydrate reserves in adult trees of
two contrasting broad-leaved species (Quercus petraea andFagus sylvatica). New Phytologist 2003, 157:605-615.
Watanabe T, Yokozawa M, Emori S, Takata K, Sumida A, Hara T:
Developing a multilayered integrated numerical model ofsurface physics-growing plants interaction (MINoSGI). Global
Change Biology 2004, 10:963-982.
Hara T: A stochastic model and the moment dynamics of thegrowth and size distribution in plant populations. J Theor Biol
1984, 110:223-239.
Hara T: Dynamics of stand structure in plant monocultures. J
Theo Biol 1984, 110:223-239.
Pacala SW, Canham CD, Saponara J: Forest models defined byfield measurements: estimation, error analysis and dynamics. Ecol Monogr 1996, 66:1-43.
Zhou G, Wang Y, Jiang Y, Yang Z: Estimating biomass and net
primary production from forest inventory data: a case studyof China's Larix forests. For Ecol Man 2002, 169:149-157.
Bonan G: Ecological Climatology. Cambridge University Press,
Cambridge; 2002.
Kurachi N: Production structure of a Japanese larch stand. In
Ph.D thesis Nagoya University; 1989. (in Japanese)
Körner Ch: Plant CO2 responses: anissue of definition, timeand resource supply. New phytologist 2006, 172:393-411.
Körner Ch, Asshoff R, Bignucolo O, Hättenschwiler S, Keel SG,
Pelaez-Riedl S, Pepin S, Siegwolf RTW, Zotz G: Carbon flux andgrowth in mature deciduous forest trees exposed to elevated CO2. Science 2005, 309:1360-1362.
Norby RJ, Ledford J, Reilly CD, Miller NE, O'Neill EG: Fine-root
procution dominates response of a deciduous forest to
atmopheric CO2 enrichment. Proceedings of the national Academyof Sciences 2004, 101:9689-9693.
Schafer KVR, Oren R, Ellsworth DS, Lai CT, Herrick JD, Finzi AC,
Richter DD, Katul GG: Exposure to an enriched CO2 atmosphere alters carbon assilination and allocation in a pine forest ecosystem. Global Change Biol 9:1378-1400.
Steffen W: The anthropocene, global change and sleeping
giants: where on Earth are we going? Carbon Balance and Manage 2006. doi:10.1186/1750-0680-1-3
Grace PR, Post WM, Hennessy K: The potential impact of climate change on Australia's soil organic carbon resources.
Carbon Balance and Manage 2006. doi:10.1186/1750-0680-1-14
Garcia-Gonzalo J, Peltola H, Briceño-Elizondo E, Kellömaki S:
Changed thinning regimes may increase carbon stock underclimate change: A case study from a Finnish boreal forest.
Climatic Change 2007, 81:431-454.
Katsuno M: A basic study about gross production in a Cryptomeria japonica plantation. In Ph.D thesis Nagoya University,
Nagoya; 1990. (in Japanese)

Page 17 of 18(page number not for citation purposes)

Carbon Balance and Management 2007, 2:6

43.
44.

45.

46.
47.
48.
49.
50.

51.
52.
53.
54.
55.

56.
57.
58.

59.
60.
61.
62.
63.
64.
65.

66.
67.

Takata K, Emori S, Watanabe T: Development of the minimal
advanced treatments of surface interaction and runoff. Global
Planet Change 2003, 38:209-222.
Sellers PJ, Randall DA, Collatz GJ, Berry JA, Field' CB, Dazlich DA,
Zhang C, Collelo GD, Bounoua L: A revised land surface parameterization (SiB2) for atmospheric GCMs. Part I: Model formulation. J Clim 1996, 9:676-705.
Barbaroux C, Breda N: Contrasting distribution and seasonal
dynamics of carbohydrate reserves in stem wood of adultring-porous sessile oak and diffuse-porous beech trees. Treephysiol 2002, 22:1201-1210.
Yokozawa M, Hara T: A canopy photosynthesis model for the
dynamics of size structure and self-thinning in plant populations. Ann Bot 1992, 70:305-316.
Shinozaki K, Yoda K, Hozumi K: A quantitative analysis of plantform-the pipe model theory, I. Basic analyses. Jap J Ecol 1964,
14:97-105.
Yokozawa M, Kubota Y, Hara T: Crown architecture and species
coexistence in plant communities. Ann Bot 1996, 78:437-447.
Lindner M, Sievänen R, Pretzsch H: Improving the simulation ofstand structure in a forest gap model. Forest Ecol Man 1997,
95:183-195.
Lincosalo T, Carter TR, Hakkinen R, Hari P: Predicting spring phenology and frost damage risk of Betula spp. Under climatic
warming: a comparison of two models. Tree Physiol 2000,
20:1175-1182.
Sarvas R: Investigations on the annual cycle of development offorest trees. Active period. Communicationes Linstinstituti Forestalis
Fenniae 1972, 76:1-110.
Sarvas R: Investigations on the annual cycle of development offorest trees. Autumn dormancy and winter dormancy. Communicationes Linstinstituti Forestalis Fenniae 1974, 84:1-101.
Tadaki Y: Leaf opening and falling of Japanese larch at different altitudes. Jpn J Ecol 1994, 44:305-314.
Cannell MGR, Smith RI: Thermal time, chill days and predictionof bud burst in Picea Sitchensis. J Appl Ecol 1983, 20:951-963.
Watanabe T, Nakai Y, Kitamura K, Utsugi H, Tobita H, Ishizuka S:
Scaling energy and CO2 fluxes from leaf to canopy using Multilayered Implementation for Natural Canopy-Environment
Relations (MINCER). Phyton Special issue: "APGC 2004" 2005,
45:353-360.
Larcher W: Physiological Plant Ecology, Ecophysiology andStress Physiology of Functional Groups. Springer; 1995.
Hozumi K, Kurachi N: Estimation of seasonal changes in translocation rates in leaves of a Japanese larch stand. Bot MagTokyo 1991, 104:25-36.
Funai R: Cambial Reactivation in Deciduous Hardwoods. In
Tree Sap Proceedings of the 2nd International Symposium on Sap Utilization (ISSU) in Bifuka 2000 Bifuka Edited by: Terazawa M. Hokkaido University Press; 2000.
Ögren E: Maintenance respiration correlates with sugar but
not nitrogen concentration in dormant plants. Physiologia
Plantarum 2000, 108:295-299.
Kujiraoka K: Seasonal changes in allometry and carbohydrate
storage of Betula ermanii in Hokkaido. In Master thesis
Hokkaido University; 2004. (in Japanese)
Kozlowski TT, Pallardy SG: Physiology of woody plants. Academic Press; 1997.
Kurachi N, Hagihara A, Hozumi K: Canopy photosynthetic production in a Japanese larch forest. II. Estimation of the canopy photosynthetic production. Ecol Res 1993, 8:349-361.
Kurachi N, Hagihara A, Hozumi K: Evaluation of the light interception by non-photosynthetic organs in a Larix leptolepis
plantation. Ecol Res 1986, 1:173-183.
Farquhar GD, von Gaemmerer S, Berry JA: A biochemical modelof photosynthetic CO2 assimilation in leaves of C3 species.
Planta 1980, 149:78-90.
Ball JT, Woodrow IE, Berry JA: A model predicting stomatal
conductance and its contribution to the control of photosynthesis under different environmental conditions. Progress in
Photosynthesis Research 1987, 4:5221-5224.
Hukami M: Measurement of root respiration of Larix kaempferi. In Graduation thesis Nagoya University; 1981. (in Japanese)
Hagihara A, Hozumi K: Estimation of canopy respiration and its
seasonal change in a Cbamaecyparis obtusa plantation. J Jpn
For Soc 1977, 59:405-413.

http://www.cbmjournal.com/content/2/1/6

Publish with Bio Med Central and every
scientist can read your work free of charge"BioMed Central will be the most significant development for
disseminating the results of biomedical researc h in our lifetime."
Sir Paul Nurse, Cancer Research UKYour research papers will be:available free of charge to the entire biomedical community
peer reviewed and published immediately upon acceptancecited in PubMed and archived on PubMed Central
yours — you keep the copyright

BioMedcentral

Submit your manuscript here:
http://www.biomedcentral.com/info/publishing_adv.asp

Page 18 of 18(page number not for citation purposes)