Interactive comment on “Simulating carbon and water cycles of larch forests in East Asia by the BIOME-BGC model with AsiaFlux data” by M. Ueyama et al

This manuscript represents a very nice study using a modeling tool to synthesize and understand the variability in carbon and water dynamics observed with eddy covariance at six larch forests located across northeastern Asia. Model biases are first evaluated and the model is then modified to minimize biases. I was particularly happy to see that he stand history was considered in the model simulations, and that it was demonstrated to be quite important. The use of the improved model was valuable for understanding spatial variability in controls over carbon and water exchange. The identification of the importance of spatial variability in allocation to above vs. belowground tissues


Introduction
The northern high latitude region is currently undergoing rapid and drastic warming (IPCC, 2007); the air temperatures in Eastern Siberia have risen by 7 • C in winter and 1-2 • C in annual basis during the past few decades (Dolman et al., 2008).The terrestrial ecosystems in this region have responded to the warming climate through various feedback Published by Copernicus Publications on behalf of the European Geosciences Union.
processes (Chapin et al., 2005;Hinzman et al., 2005).These changes are likely to affect the carbon and water cycles in boreal and cool temperate forests, altering snow cover, permafrost dynamics, growing season length, and the severity of summer drought.
The predicted climatic changes affect the seasonal weather patterns differently (e.g.Manabe et al., 1992), and it is particularly important to understand how the changes in weather patterns affect terrestrial carbon and water cycles.The effects of current and future climatic changes on the terrestrial carbon and water cycles over northern high-latitude ecosystems are complex.Detected and/or potential effects include: (1) the earlier onset of photosynthetic activity in high-latitude forests due to the spring warming (Myneni et al., 2001;Kimball et al., 2004) leads to greater productivity and thus potentially increases the carbon sink (e.g.Keyser et al., 2000;Euskirchen et al., 2006;Welp et al., 2007); (2) the warming climate potentially stimulates the decomposition of soil carbon (Valentini et al., 2000;Piao et al., 2007;Ueyama et al., 2009) and reduces the carbon sink; and (3) atmospheric drying over the regions of warming could affect the carbon and water cycles of boreal forests since the carbon and water fluxes of boreal forests in Eurasia are highly sensitive to changes in atmospheric humidity (e.g.Schulze et al., 1999;Wang et al., 2005;Ohta et al., 2008).
A number of terrestrial biosphere models have been applied to simulate the carbon and water cycles in high-latitude ecosystems (Euskirchen et al., 2006;Kimball et al., 2007).Although these models were tested with the data observed over the high-latitude biomes (Amthor et al., 2001;Grant et al., 2005), these validations were mostly conducted for evergreen conifer (Cienciala et al., 1998;Clein et al., 2002;Ueyama et al., 2009), deciduous broadleaf (Kimball et al., 1997a), or arctic tundra ecosystems (Engstrom et al., 2006), but rarely conducted for the deciduous conifer (larch) forests.The lack of validation studies for larch forests could induce uncertainties in predicting the high-latitude and global carbon and water cycles.
Larch forests are widely distributed over many cooltemperate and boreal regions in the northern hemisphere (Gower and Richards, 1990).Across Eurasia to East Asia, a number of eddy covariance measurements have been conducted for several larch forests near the arctic (e.g.Ohta et al., 2001;Machimura et al., 2005;Nakai et al., 2008) and over boreal (Wang et al., 2005;Li et al., 2005) and cool temperate (Hirata et al., 2007) regions.These measurements have revealed the important processes that determine the carbon and water cycles in larch forests, such as the environmental factors that control evapotranspiration (Ohta et al., 2008) and carbon flux (Hollinger et al., 1998;Wang et al., 2005;Li et al., 2005;Hirata et al., 2007;Nakai et al., 2008), and the role of stand disturbance (Machimura et al., 2007).Since these analyses are site-specific, we need to analyze them at multiple site scales to clarify how the carbon and water fluxes are spatially distributed, how the responses to the environmental conditions differ spatially, and how terrestrial biosphere models perform in simulating carbon and water dynamics.These studies are limited for larch forest ecosystems because of few observations in comparison, compared with the other forest ecosystems in North America (Baldocchi et al., 2001;Thornton et al., 2002) and Europe (Valentini et al., 2000).However, the recent availability of observations makes it possible to conduct synthesis studies for this ecosystem.
This study is a first step to understand the carbon and water cycles of larch forest in northern Eurasia to East Asia.In this study, a process-based terrestrial biosphere model, BIOME-BGC (Thornton et al., 2002), was tested for larch forests at six AsiaFlux sites and used to identify important environmental factors that affect the carbon and water cycles temporally and spatially.Our specific objectives are to (1) improve the model performance by using the observed flux data; (2) clarify the environmental factors, including climate and stand disturbance, that control the carbon and water fluxes across the larch forests in East Asia; and (3) examine the response of the carbon budget to changes in seasonal weather patterns.

Site descriptions
The present study is a synthesis of 12 years of site data from six larch forests from Eurasia to East Asia (Fig. 1).The sites include the Tomakomai site (TMK), Japan; the Laoshan site (LSH), China; the Southern Khentii Taiga site (SKT), Mongolia; the Yakutsk site (YLF), Russia; the Neleger site (NEL), Russia; and the Tura site (TUR) in Russia.These sites are part of the AsiaFlux network (Hirata et al., 2008;Saigusa et al., 2008;Mizoguchi et al., 2009), covering a broad range of climate with annual precipitation totals from 240 mm to 1750 mm per year and annual mean air temperatures from −10 • C to 10 • C. In terms of seasonal climatic variations (Fig. 2), TMK is characterized as a cool temperate climate with a cool early summer, warm winter and spring, and high humidity over the course of the year (low summer VPD and sufficient precipitation).The LSH and SKT sites are characterized by a boreal climate with cold winters, hot summers, dry atmospheric conditions, although LSH has much more precipitation during the summer than SKT.All of the Russian sites (YLF, NEL, and TUR) are characterized by a severe continental climate with cold winters, hot and dry summers, and low precipitation totals.Details of site information are summarized in Table 1.
TMK is a planted Japanese larch (Larix kaempferi) forest located at approximately 15 km northwest of Tomakomai, Hokkaido, in northern Japan (Hirano et al., 2003;Hirata et al., 2007).The stand was planted in 1957-1959, and the stand age was about 45 years old at the time of this study.The forest was sparsely dominated by broadleaf    trees (Betula ermanii, Betula platyphylla, and Uhnus japonica) and spruce (Picea jezoensis) with an understory of ferns (Dryopteris crassirhizoma and Dryopteris austriaca).The maximum LAI of larch, other overstory species, and understory species were about 3.1, 2.5, and 3.6 m 2 m −2 in summer, respectively.
LSH is a planted larch (Larix gmelinii) forest in northern China (Wang et al., 2005a, b).The forest was established as a plantation in 1969 as a pure larch forest.The maximum LAI of larch trees was about 2.5 m 2 m −2 .Some broadleaf species sparsely consisted of the canopy (e.g.Betula platyphylla, and Fraxinus mandshurica).The forest floor was covered by the understory vegetations, shrubs, and species in the Cyperaceous and Liliaceous families.
SKT is a Siberian larch (Larix sibirica) forest located at approximately 25 km northeast of the Mongonmorit village in the Tov province of Mongolia (Li et al., 2005(Li et al., , 2006(Li et al., , 2007)).The age structure of the larch ranged from ∼ 70 to over 150 years old, but the oldest trees were older than 300 years.The forest experienced a large-scale fire in 1996-1997 in which approximately 37.5% of the total stem density was burnt.The understory is composed of a distinct layer of grass and scattered shrubs.The maximum LAI of the larch and understory species were about 2.7 and 1.7 m 2 m −2 , respectively.
YLF is a Cajander larch (Larix cajanderi) forest located approximately 20 km north of Yakutsk, Russia, eastern Siberia (Ohta et al., 2001(Ohta et al., , 2008;;Dolman et al., 2004).The larch trees are on permafrost soil with an active layer depth of approximately 1.2 m.The LAI of the larch trees was 1.56 m 2 m −2 , whereas the forest floor was covered by a dense understory of vegetation, such as Vaccinium species, with an LAI of 2.1 m 2 m −2 .Since there has not been a fire at the site for at least the last 80 years, the average stand age is 160 years.www.biogeosciences.net/7/959/2010/Biogeosciences, 7, 959-977, 2010 NEL is a larch (Larix cajanderi) forest located on the approximately 30 km northwest of Yakutsk, Russia (Machimura et al., 2005(Machimura et al., , 2007)).The forest is also on continuous permafrost with an active layer depth of approximately 1.0 m.The average tree height was 8.6 m, and the maximum tree height was approximately 21 m.The forest floor was covered by Vaccinium uliginosum and Pyrola incarnate with moss species.The LAI of the larch trees was 1.9 m 2 m −2 , whereas that of the understory was 1.0 m 2 m −2 .
TUR is a Gmelin larch (Larix gmelinii) forest located in Tura in the Evenkia Autonomous District in Central Siberia (Nakai et al., 2008).The stand age is about 105 years old.The LAI of the larch trees was ∼ 0.3 m 2 m −2 .The soil type is Cryosol, and the permafrost table is within 1 m depth.The forest floor was covered by some shrub species, such as Betula nana, Ledum palustre, Vaccinium uliginosum, and Vaccinium vitis-idaea.The ground was densely covered with lichens and mosses.

Data processing
We used the observed carbon and water fluxes from the eddy covariance method for the six tower-sites.The observed net ecosystem exchange (NEE) were partitioned into ecosystem respiration (RE) and gross primary productivity (GPP) using relationships obtained from the nighttime NEE and air temperature (Hirata et al., 2008;Takagi et al., personal communications, 2009).Quality control and gap filling of the data were conducted by standardized methods (Hirata et al., 2008;Takagi et al., personal communications, 2009).Observed half-hourly fluxes and meteorology were aggregated to a daily basis for the model input.In this study, negative and positive NEE values indicate a net sink and source of atmospheric carbon, respectively.

BIOME-BGC model
The carbon and water cycles of the larch forests were simulated using the BIOME-BGC model (Thornton et al., 2002).The BIOME-BGC is a process-based terrestrial ecosystem model that simulates biogeochemical and hydrologic processes across multiple biomes.The model is driven at a daily time scale with the meteorological values of the daily maximum, minimum, and average air temperatures, precipitation, daytime VPD, and solar radiation.
Evapotranspiration (ET) is estimated from the sum of transpiration, evaporation from the soil and canopy, and sublimation (Kimball et al., 1997b).Evapotranspiration is basically calculated by a Penman-Monteith model (Running and Coughlan, 1988).Soil conductance is simply estimated from the number of days since a rainfall event.Canopy evaporation is estimated from the amount of water intercepted by the canopy.Transpiration is regulated by the canopy conductance under the daily meteorological conditions of the minimum air temperature, VPD, solar radiation, and soil water availability (Jarvis, 1976).
The model has three compartments for carbon and nitrogen: vegetation, litter, and soil.Each compartment is subdivided into four pools on the basis of differences in their function (e.g., leaf, stems, coarse roots, and fine roots) and residence time (e.g., active, intermediate, slow, and passive recycling).GPP is estimated by coupling the Farquhar biochemical model (Farquhar et al., 1980) with a stomatal conductance model (Jarvis, 1976).Ecosystem respiration (RE) is calculated as the sum of autotrophic and heterotrophic respiration (AR and HR, respectively).The AR and HR are calculated from the carbon and nitrogen pools and the temperature (for AR and HR) and soil water condition (for HR only).Further details for the BIOME-BGC model have been described in previous papers (e.g., Kimball et al., 1997a, b;Thornton et al., 2002;Ueyama et al., 2009).

Model modifications
To improve the model performance for the larch forests, we applied three modifications to the original model.The requirement of these modifications were found by checking the default model output with field observations of carbon and water dynamics, seasonal variation in snow cover, and soil temperature.First, we incorporated the seasonality of the percent of leaf nitrogen in Rubisco (PLNR) based on the strong seasonal variation in the leaf C:N ratio in larch trees (Kim et al., 2005) and the seasonality of PLNR in deciduous trees (Wilson et al., 2000) as follows: where PLNR base is the maximum PLNR (which corresponds to the PLNR of the original BIOME-BGC model); CN leaf is the leaf C:N ratio; CN litter is the litter C:N ratio; d is the day of the year; d onset is the onset day; and d growth is the length of the growing period.The onset day and growing season length are calculated by the BIOME-BGC model based on an empirical phenology model (White et al., 1997).The incorporation of this seasonality produces a higher GPP in the beginning of the growing season and then a gradual decline of GPP later in the growing season, results which are consistent with flux tower-based GPP measurements (Hirata et al., 2007).Second, the original BIOME-BGC model substantially overestimated the snow sublimation, which induced an unreasonable snow disappearance in late winter and thus a water deficit in the boreal summer and autumn.We simply modeled the snow disappearance by using only the daily air  Table 2 a Determined each site.b Kajimoto et al. (1999) c 0.045 for the sites of TMK, SKT, YLF, NEL, and TUR, whereas 0.005 for LSH.d Matyssek and Schulze (1987) e ranged in previous examined values: Pietsch et al., (1991), andVygodskaya et al. (1997) f Turner et al., (2003) g Pietsch et al. (1991) Parameters used in the BIOME-BGC model.* PSI, leaf and soil water potential * SLA, specific leaf area a Determined each site.b Kajimoto et al. (1999) c 0.045 for the sites of TMK, SKT, YLF, NEL, and TUR, whereas 0.005 for LSH.d Matyssek and Schulze (1987) e ranged in previous examined values: Pietsch et al. (1991), andVygodskaya et al. (1997) f Turner et al. (2003) g Pietsch et al. (1991) temperature in order to make the model reproduce the observed variations of water and carbon fluxes.
Third, we used the monthly moving average for the estimation of soil temperature for further reproduction of the observed results, although the original BIOME-BGC model estimated the soil temperature by using the 10-day moving average of the daily air temperature.The estimated soil tem-perature during the snow period was also overestimated in the TMK site by the original BIOME-BGC model; this was probably caused by the overestimation of the snow insulation effect probably because the difference in snow depth, snow density, and soil porosity among the sites.We changed the parameter for the insulation effect from 0.83 to 0.55 at the simulations of TMK to reproduce the observed soil www.biogeosciences.net/7/1/2010/Biogeosciences, 7, 1-19, 2010 temperature.This calibration decreased the soil temperature during the snow period and thus decreased the root respiration and HR in the winter.

Model initializations
The model was initialized by a spinup run, in which the dynamic equilibriums of the soil organic matter and vegetation components were determined by using a constant CO 2 concentration of 296 ppm and the observed daily meteorological data and eco-physiological parameters (Table 2).We used soil texture data (Zobler, 1986) and rooting depth data (Webb et al., 1993) as the model inputs for the sites where site-specific observation is not available.Only for the YLF, NEL and TUR sites, we used the observed seasonal maximum of the active layer depth as effective rooting depth.Using the spinup-endpoint as an initial condition, we conducted the model simulations from 1915 to the final year of the observations with the ambient CO 2 concentration (Enting et al., 1994;Tans and Conway, 2005).In these simulations, we repeatedly used the observed meteorology at the site rather than the long-term meteorological records from the weather stations near the sites in order to avoid biases due to data discontinuities.Since the observations at YLF, NEL, and TUR were only conducted during periods of the vegetation growth, we used the meteorological data of the National Climate Data Center (NCDC; Global Surface Summary of Day version 8) from weather stations near the sites for the model input for those winter periods.We implemented simple disturbance treatments for TMK and LSH as plantation and SKT as fire based on the definitions presented by Thornton et al. (2002).For the plantation disturbance of TMK and LSH, we assumed that the amount of vegetation planted was 30% of the dynamic equilibrium of the vegetation, and the rest was removed from the sites; the leaf and fine root C and N pools were included in the fine litter pool.The aboveground live and dead wood C and N pools were removed from the site, and the belowground live and dead wood C and N pools were included in the coarse woody debris (CWD) pools (Thornton et al., 2002).Based on the stand history (Hirano et al., 2003;Wang et al., 2005a), we assumed that the plantations were established in 1959 for TMK and in 1969 for LSH.Since the stands experienced a large-scale fire in 1996 and 1997 at SKT (Li et al., 2005), we also implemented a fire disturbance in which 37.5% of the forest was affected.The proportions of the leaf, fine root, live wood, and fine litter C and N pools were assumed to be consumed in the fire and lost to the atmosphere, and the affected portions of the dead wood C and N pools were sent to the CWD pools (Thornton et al., 2002).Although the stand disturbance history of the mature forest sites in YLF, NEL, and TUR were not available, we inserted a fire disturbance just after the spinup run (year of 1915) for YLF, NEL, and TUR, assuming that the same amount as in SKT, 37.5%, was affected by the ground fire based on the fact that Siberian larch forests experience a ground fire every 20-50 years (Schulze et al., 1999).

Model applications
We first simulated the carbon and water fluxes by using the original BIOME-BGC model with the original parameters for deciduous needle leaf forest (DNF) shown in Table 2 (White et al., 2000).Then, we calibrated the modified Biogeosciences , 7, 959-977,      model by using the AsiaFlux data observed at the six larch forests, and evaluated how the AsiaFlux data could improve the model performance for the larch forests in East Asia.
Most of the default ecophysiological parameters for DNF are from those for evergreen needleleaf forest (White et al., 2000) owing to limitation of data availability.We updated the new general parameters for DNF in Table 2, referring the literature available data; all of parameters were common across the six sites except for the two site-specific parameters.The model was parameterized according to the following criteria: (1) the parameters were at first changed at a minimum in order to reduce biases in the fluxes, such as the magnitude of day-to-day variations and unrealistic declines due to suppression effects, by referring to the previous studies; and (2) the peak of GPP was adjusted by changing the allocation ratio between new fine roots C to new leaf C (A root_leaf ).
Assuming that the allocation ratio between aboveground and belowground components strongly depends on the climatic conditions (e.g., Friedlingstein et al., 1999;Chapin et al., 2002;Vogel et al., 2008), the parameter for A root_leaf was determined for each site.The canopy interception coefficient in LSH was only changed to 0.005 LAI −1 d −1 in order to reproduce the observed ET.In both the original and the improved simulation, we applied the stand disturbance treatment in the model initialization (as described in Sect.3.3).To understand the role of stand disturbance on the annual carbon budget, the observed annual carbon budgets were also compared to the simulated results with and without the stand disturbances.
Next, the spatial distributions of the modeled carbon fluxes were examined by using the annual climatology of each site.The simulated annual values of carbon fluxes from the improved model were compared with the annual climatology of air temperature, precipitation, and radiation to clarify the important climatic variables for determining the spatial variations in the fluxes.
Then, to understand the response of the carbon cycle to changes in the weather conditions, we conducted a sensitivity analysis with biased meteorological data, by giving either warmer or cooler, wetter or drier, and sunnier or cloudier biases for each season; winter (December to February), spring (March to May), summer (June to August), and autumn (September to November).The biases were determined as 3σ (σ is the standard deviation) of the long-term natural climate variabilities (Table 3).The natural climate variabilities for air temperature, precipitation, and VPD in TMK were estimated by using the weather station data observed by the Japan Meteorological Agency (JMA) from 1961 to 2000, whereas those in LSH, SKT, YLF, NEL, and TUR were estimated from the observations from NCDC from 1981 to 2000.The variabilities for radiation were derived from the International Satellite Cloud Climatology Project Radiation Flux Profile data from 1984 to 2000 (ISCCP-FD; Zhang et al., 2004) centered on the observed site.

Model validation
The original model generally underestimated carbon flux variations at daily, monthly, and annual time scales (Figs. 3,5 and 6).At the daily time scale, the results were not consistent with the observed fluxes.They showed considerable scatter and underestimated or overestimated the observed values in each site (Fig. 3).The slopes between the observations and the simulation results were 0.65 for NEE, 0.73 for GPP, 0.60 for RE, and 0.70 for ET, which all deviated substantially from 1.0, and those of R 2 were 0.53, 0.75, 0.72, and 0.25, respectively.The monthly variations of the simulated and observed fluxes also showed that the amplitude of the flux variation simulated by the original model was inconsistent with the observed results (Fig. 5).The peaks of GPP and RE were generally underestimated for TMK, SKT, and YLF and overestimated for NEL and TUR, although the seasonality, such as the onset and offset of the growing season, was generally reproduced.The annual GPP, RE, and ET were underestimated in TMK, but overestimated in other sites, resulting in those slopes of regression lines that deviated from 1.0 (Fig. 6).
Using the observed fluxes as calibration data, the improved model successfully simulated the carbon and water fluxes at daily, monthly, and annual time scales (Figs. 4,5 and 6).The seasonality of the carbon fluxes was also substantially improved (Fig. 5).Although the observed peak of GPP and NEE in the TMK site in June could not be reproduced by the original model, the improved model successfully reproduced the phase and amplitude of the carbon fluxes owing to the incorporation of the seasonality of PLNR, which was consistent with the observed results at the leaf (Kim et al., 2005) and canopy scales (Hirata et al., 2007).This improvement suggested that using a constant nitrogen ratio in the leaf could contribute to the uncertainty in the original BIOME-BGC model.In the boreal sites of YLF, NEL, and TUR, since the snow sublimation was reduced due to the modification, the water conserved during the winter was used in the summer periods in the improved model and thus mitigated the water stress (Fig. 5d), resulting in the improvement in an seasonality of ET for YLF.Although the carbon fluxes were substantially improved in the model, the simulated ET at the daily time scale was not substantially improved, compared with the carbon fluxes.This discrepancy might be caused by both the observation and simulation (shown in Sect.4.4).
The improvements of the model reduced the root mean square error (RMSE) at almost every site (Table 4), and the slops became closer to 1.0 (Figs. 4 and 6).The results from the improved model simulation highly correlated with the observed GPP (R 2 = 0.99), RE (R 2 = 0.98), and ET (R 2 = 0.86), compared to the results from the original simulation (Fig. 6).The annual carbon budget simulated by the model was also improved by the calibrations.On the other  hand, the improved model still only explained 41% of the variation in the observed NEE carbon budget (Fig. 6a); the simulated results by the improved model had a tendency to underestimate the observed carbon uptake.
The inclusion of the stand history was an important factor for improving annual NEE estimation (Fig. 7).If we run the model without including stand history, the simulated NEE was substantially underestimated in both the original and improved models.The simulations run without including the stand history could not reproduce the observed annual NEE (R 2 = 0.13), whereas those considering the stand history reproduced 50% of the variation.The model performance for estimation of the annual NEE was improved by considering the stand history; the annual sink was increased in the planted sites of TMK and LSH and the fire-disturbed site of SKT by considering the stand histories.
Overall, the model was successfully improved by the calibration with regional network of flux tower measurements.The improved model reasonably reproduced the carbon and water fluxes at the daily, monthly, and annual time scales, although the daily water flux still had uncertainties.We used this improved model to evaluate the spatial distributions of the fluxes and the sensitivity of the fluxes to the weather conditions.

Spatial gradient of carbon cycle and climate
The sensitivities of the simulated annual carbon and water fluxes to the annual climate are shown in Fig. 8. Across the six sites, the simulated GPP and RE substantially increased with the annual air temperature (R 2 = 0.84; p < 0.01 for GPP and R 2 = 0.82; p < 0.01 for RE) as well as the annual precipitation (R 2 = 0.86; p < 0.01 for GPP and R 2 = 0.85; p < 0.01 for RE), and there was no clear relationship between those fluxes and the amount of radiation (R 2 = 0.20; p = 0.38 for GPP and R )

Figure 8
Relationships between mean annual climate variables (air temperature, precipitation, and radiation) and simulated annual GPP, RE, and NEE.Symbols are as in Table 1.  1.
there was a strong correlation between annual air temperature and precipitation across the sites, we could not identify which of these processes was the more important controlling factor for explaining the spatial variation of those fluxes.Since the simulated annual fluxes were consistent with the observed results (Fig. 6), these trends were also applicable to the observed fluxes.Our results indicated that the air temperature and/or precipitation are potentially important controlling factors for explaining the spatial distribution of the carbon and water fluxes of larch forests.The simulated NEE was weakly correlated with the annual climate (Fig. 8) in that the annual carbon sink increased with an increase of the annual mean air temperature (R 2 = 0.55; p = 0.09), precipitation (R 2 = 0.46; p = 0.14), and radiation (R 2 = 0.50; p = 0.12).Based on the correlation coefficients between NEE and the climate variables, the annual mean air temperature was most strongly related to the spatial distribution of NEE.Although both the GPP and RE had positive correlations with the climate variables, a greater GPP induced a greater annual carbon sink, whereas a greater RE did not reduce this sink.This result suggests that the spatial distribution of NEE was generally determined by that of GPP rather than RE.Since the annual carbon budget was highly sensitive to the stand history (Fig. 7), the spatial distribution of NEE was determined by both climate variables and stand disturbance history.
In the model calibration, we adjusted the seasonality of the fluxes by changing the parameter for the allocation ratio of carbon between aboveground and belowground (A root_leaf ).The estimated site-specific A root_leaf showed a weak negative correlation with the mean annual air temperature (R 2 = 0.24; p = 0.30) and total precipitation (R 2 = 0.24; p = 0.18) (Fig. 9).These trends became statistically significant (R 2 = 0.91; p = 0.02 for air temperature, R 2 = 0.77; p = 0.05 for total precipitation) when the relationships were estimated without LSH (dashed lines in Fig. 9), since the carbon fluxes of LSH might be affected by tree species other than larch due to the limited footprint (H.Wang, personal communication, 2009).The negative trends indicate that the forests could allocate more carbon aboveground under warmer and wetter climate.The estimated trend of A root_leaf can be explained by the allocation strategy of plants that adjust resource acquisition to maximize the capture of the most limiting resource (Chapin et al., 2002); light could be the most limiting resource under a favorable climate, whereas soil water and other belowground resources might limit net primary productivity (NPP) under a cold and dry environment.A similar re-sult was reported by Vogel et al. (2008), who summarized the carbon allocation in boreal black spruce forests across North America; soil temperature was positively correlated with the aboveground NPP and negatively correlated with the belowground NPP across those sites.The higher allocation to belowground at northern sites is also consistent with the field survey of root systems at TUR. Kajimoto et al. (2003) observed the greater rooting area of the larch trees compared with the crown projected area, and concluded that the larch trees on permafrost soil competed for belowground resources rather than light.The simulated LAI was negatively correlated with the A root_leaf (R 2 = 0.83) (Fig. 9), suggesting that A root_leaf strongly constrained the LAI in the model simulation.The climatic condition apparently strongly controls the allocation strategy (Fig. 9a, b) and thus controls the carbon fluxes (Fig. 8) through effects of allocation on LAI (Fig. 9c).The strong correlation between A root_leaf and the simulated LAI indicates that the satellite-derived LAI might be useful to infer the spatial distributions of A root_leaf and improve the simulation of regional carbon dynamics.This allocation parameter might also be inferred from the climate data used to drive the regional simulations.Such a dynamic allocation scheme has been incorporated in some biosphere models (e.g.Friedlingstein et al. 1999;Sitch et al., 2003) and will likely improve the performance of the future projections made with the BIOME-BGC model.

Effect of seasonal climate anomalies on the carbon cycle
To understand the responses of carbon cycle to the weather conditions of air temperature, precipitation, VPD, and solar radiation, we examined the sensitivity of the larch forests to changes in the seasonal weather patterns (Sects.3-4).The results of the sensitivity analysis are shown in Table 5. Spring warming enhanced the carbon sink (Table 5 case ta), but the summer and autumn warming tended to decrease the sink www.biogeosciences.net/7/959/2010/Biogeosciences, 7, 959-977, 2010 (Table 5 case tb).The spring warming enhanced GPP rather than RE except at SKT, where it resulted in the increase of the annual carbon sink (Table 5 case ta).This was because the spring warming induced the earlier onset of leaf emergence and enhanced photosynthetic activity, whereas RE was well regulated by the low soil temperature.On the other hand, summer, autumn, or winter warming decreased the annual carbon sink, except in the autumn of NEL, because the warming stimulated RE more than GPP (Table 5 case tc, te and tg).These results were consistent with the tower observations at TMK (Hirata et al., 2007), SKT (Li et al., 2005), and another Siberian site (Hollinger et al., 1998).Tree-ring data for the Siberian larch also showed that these larch trees did not benefit from mid-summer to autumn warming (Vaganov et al., 1999;Kirdyanov et al., 2003).These results suggest that the future carbon budget of larch forests depends on the changes in warming during the spring, summer, and autumn periods.
A change in precipitation substantially affected GPP and RE in SKT (Table 5 case pc and pd) in that a wetter summer enhanced both GPP and RE by mitigating the water stress to photosynthesis and decomposition.For other sites, a change in precipitation did not substantially influence GPP and RE, but it greatly influenced NEE, especially for the boreal sites of YLF, NEL, and TUR.The high sensitivity of NEE occurs because NEE is a small difference between two large fluxes.The forests of TMK and LSH did not respond strongly to the change in precipitation (Table 5 case pa-ph), a result which was consistent with the observations (Hirata et al., 2007).
The sensitivity of the annual carbon budget to the change in VPD was only obvious in the boreal summer (Table 5 case vc and vd).The higher VPD in the summer regulated GPP and thus decreased the annual sinks of the LSH, SKT, YLF, and TUR sites (Table 5 case vc).On the other hand, the lower VPD in the summer mitigated the regulation of photosynthesis at SKT, resulting in the increase of the annual sink there (Table 5 case vd).This sensitivity to VPD was consistent with the observed results for LSH (Wang et al., 2005b), SKT (Li et al., 2005), and TUR (Nakai et al., 2008) as well as a comparative study between the TMK and LSH sites (Wang et al., 2005b).The reason for the low sensitivity to VPD at the TMK site is that the summer VPD is not high enough to sufficiently reduce the stomatal conductance (Fig. 2).
The GPP in the TMK site was the most sensitive to the change in solar radiation.An increase of approximately 12% of the annual GPP resulted in an approximately 70% enhancement of the annual carbon sink, compared to the control simulation (Table 5 case sc).The decline of summer radiation greatly decreased GPP and the magnitude of negative NEE, meaning a decrease in the annual sink at the TMK site (Table 5 case sd).Among the boreal sites of LSH, SKT, YLF, NEL, and TUR, the change in the radiation did not substantially affect the GPP regardless of the season, but the annual carbon sink (negative NEE) was positively correlated with the change in summer radiation.This higher sensitivity to radiation was also reported for the observed results at TMK (Hirata et al., 2007).Considering the sensitivity to VPD, the summer VPD could be an important factor that controls the variations of the carbon budget, but the summer radiation might become an important factor if the regulation by VPD is not strong.
The sensitivity studies identified the important weather conditions that control the annual carbon cycle in each site.In the temperate site of TMK, summer radiation was the most important control over the simulated annual carbon budget.On the other hand, summer VPD was most important in controlling the simulated annual carbon budgets of boreal sites LSH, SKT, and YLF.Air temperature was also an important controlling factor of the annual carbon budgets among all the larch forests; a warmer spring could enhance the carbon sink, but a warmer summer and autumn could decrease the sink.This temperature sensitivity was the most important controlling factor of the simulated carbon budget for the TUR and NEL sites.The sensitivity study also indicated that the small changes in GPP and RE induced large shift in the annual carbon balance, and the balance could be quite sensitive to the seasonal weather anomalies.

Further model improvements and potential limitations
This study, which is a first attempt to validate a terrestrial biosphere model at multiple larch forests sites in northern Eurasia to East Asia, identified several biases in the default model.Our analysis demonstrated that the calibration by the AsiaFlux observations greatly improved the BIOME-BGC model at daily, monthly, and annual time scales, and the improved model enabled us to perform more accurate sensitivity studies.However, we should note the following limitations and potential improvements.
The field-estimated GPP and RE were derived from NEE measurements during nighttime period.Nighttime fluxes by the eddy covariance technique have been recognized to contain some uncertainties (e.g., Gu et al., 2005;Papale et al., 2006).In addition to the uncertainties from nighttime measurement, it is difficult to accurately partition NEE into GPP and RE because there is no consensus on the partitioning method (Reichstein et al., 2005;Richardson and Hollinger, 2005).Although this study standardized the partitioning methodology among the sites, further studies will be required to estimate GPP and RE from field observations.The simulated ET was still inconsistent with the observed ET at a daily time scale (Fig. 4).This disagreement was mainly caused by the intercepted evaporation term; the simulated intercepted evaporation was substantially larger during the rainy days, compared to the observed ET.Considering that the eddy covariance method could not perfectly measure ET during and just after rainfall (e.g.Kosugi et al., 2007), we could not validate the simulated intercepted evaporation.Ohta et al. (2001) reported that the intercepted evaporation was about 15% of the precipitation at the YLF site, a result which is consistent with the simulated intercepted evaporation value of about 18% of the precipitation.The parameter for the canopy interception was needed to be site specific for the LSH site, which was probably caused by the limited footprint in this site (H.Wang, personal communication, 2009).Observed ET by the eddy covariance measurements also has uncertainties due to the energy imbalance problem (Wilson et al., 2002).The energy balance ratio at our study sites ranged from 66% to 100% (Hirata et al., 2005;Liu et al., 2005;Wang et al., 2005;Nakai et al., 2008;Ohta et al., 2008), suggesting significant uncertainties in observed ET and potentially in the field-estimated carbon fluxes.Consequently, further investigation for both observation and modeling might be required to reveal the cause of the inconsistency.
Although the simulated sensitivities to air temperature, VPD and radiation were consistent with the observed studies, the simulated response to precipitation still had large uncertainties.At the sites of TMK and LSH, the amount of precipitation did not substantially affect the modeled carbon budget; which was consistent with observations.On the other hand, precipitation substantially affected the carbon budgets in the northern sites of SKT, YLF, NEL, and TUR.Although the discrepancy in NEE might be exaggerated by small errors of the two large terms, the errors were partly caused by a lack of important processes in the model, such as the water supply from the deeper layer that the larch trees use when the water supply in the upper soil layer is limited (Li et al., 2006).Studies at YLF (Ohta et al., 2008) and NEL (Machimura et al., 2007) showed a 1-year lag of the carbon and water fluxes to precipitation probably due to the availability of water stored in the permafrost.Consequently, further improvements of the model may be required to properly represent the response to water availability.Other factors that might contribute to the discrepancy may be understory vegetation and/or permafrost dynamics (Ueyama et al., 2009).
The BIOME-BGC model considers only tree species in ecosystem, but larch stands generally consist of sparse canopy, coexisting with other species.Our study sites also contain understory vegetation other than larch (shown in Sect.2.1), which may influence carbon dynamics.
Our simulation showed that the stand history was the most important factor for predicting the carbon budget (Fig. 7).www.biogeosciences.net/7/959/2010/Biogeosciences, 7, 959-977, 2010 Although the spatial distribution of NEE was also highly sensitive to the climate variables (Fig. 8), the observed large carbon sink of TMK and SKT could not be reproduced without the stand history.These results suggest that the availability of information on stand history could improve the carbon budget simulation both at stand (the YLF, NEL, and TUR sites, where stand history was not available) and regional scales.
Several model assimilation studies point out the problem of model equifinality in that many different model parameterizations are able to reproduce the observed data (e.g.Franks et al., 1997;Schulz et al., 2001).To minimize the problem of the model equifinality, we (1) modified the submodels by carefully checking the default model output with field observations other than the carbon and water fluxes, (2) examined the high sensitive parameters across the multiple sites through the sensitivity analysis, and (3) updated the parameters mainly from the values available in the literature.In addition to validation of the current status of the fluxes, we also validated the ecosystem sensitivity to the seasonal weather anomalies, and obtained a reasonable sensitivity compared with the observed studies in each site.Consequently, the parameters are applicable to simulate the current status of larch forests of the multiple sites.While the model is well validated in this study, further data acquisition will improve the model parameterizations.
In this study, we initialized the model with the short-term meteorology at the tower sites rather than the long-term meteorological records from the weather stations near the sites, in order to avoid biases due to data discontinuities.Although the meteorological years used in this study did not substantially deviate from the long-term average (Nakai et al., 2008;Ohta et al., 2008;Hirata et al., 2007), this initialization might lead to some biases in the simulation.Since the historical changes in the climate affect the carbon and water cycles through various indirect feedbacks, such as changes in the pool size (Keyser et al., 2000;Ueyama et al., 2009), the use of short-term meteorology might ignore these effects and lead to simulations that were more indicative of the steady state condition compared with the actual carbon and water cycles.The underestimation of the annual carbon sink by the model simulation (Fig. 6) might also be caused by this initialization process.Since the simulation of sensitivity to the seasonal weather conditions (Table 5) has been conducted for a short-term change (interannual time scale), the response of the carbon and water cycles to the long-term climate anomalies (e.g.decadal or longer time scales) should be examined in future studies.

Conclusions
The regional network of tower flux measurements, AsiaFlux, helped improve the model performance of the BIOME-BGC for larch forests in East Asia.Using the flux data at six sites from Siberian subarctic to cool temperate regions, we deter-mined new general parameters for larch forests.The use of these validated parameters will allow us to conduct regionalscale simulations.The site-specific parameter of the allocation ratio of new fine root C to new leaf C was negatively correlated with the annual climatology of air temperature and precipitation; this strongly constrained the simulated LAI.In the regional scale simulation, the allocation parameter will be estimated from climate data or satellite derived LAI.
The model simulations showed that the spatial variations of GPP and RE were correlated with the variations in the annual climatology of air temperature and/or precipitation.In addition to the climatic conditions, the stand history was another important factor for explaining the spatial variations of NEE.These simulation results suggest that the spatial map of the disturbance history, such as plantation areas and fire, will be important for the projection of the regional scale carbon budget.For the future projections of the carbon cycle, the BIOME-BGC model will necessarily be coupled with a sophisticated disturbance algorithm at regional and continental scales.
The sensitivity analysis showed that spring warming could enhance the carbon sink, but summer and autumn warming could decrease the carbon sink; this is consistent with the observed results.In one temperate site, where water conditions did not regulate the carbon cycle, summer radiation was the most important factor that controlled the annual carbon budget.On the other hand, summer VPD strongly regulated the annual carbon sink in the boreal sites.Our analysis has allowed us to understand sensitivity at seasonal to interannual time scales, but it is not clear how long-term changes in climate affect the carbon and water cycles through various feedback processes.These projections should be examined in ongoing and long-term observations.
Through comparison to observations, the model uncertainties have been addressed in this study, especially for boreal regions.It is possible that the modeled response to the water conditions is inaccurate.In some boreal sites, the regulation of stomatal conductance and heterotrophic respiration by the water deficit might be overestimated.These discrepancies could be due to various causes, such as water storage in deeper soil layers, the effects of understory vegetation, and permafrost dynamics.Additional validation with available chamber and biometric measurements will help reduce the model uncertainties.

Fig. 1 .
Fig. 1.Location of study sites from northern Eurasia to East Asia.The gray area represents a distribution of larch forest based on the MOD12 landcover classification.

Fig. 2 .
Fig. 2. Seasonal variations of climate variables: (a) air temperature, (b) daytime VPD, (c) precipitation, and (d) daytime solar radiation.The data are derived from weather stations near the site from 1961 to 2000, except that the daytime solar radiation are from the International Satellite Cloud Climatology Project Radiation Flux Profile data(Zhang et al., 2004) from 1984 to 2000.The climate data of TMK are from the Japan Meteorological Agency, and those of other sites are from the National Climate Data Center.

Fig. 8 .
Fig.8.Relationships between mean annual climate variables (air temperature, precipitation, and radiation) and simulated annual GPP, RE, and NEE.Symbols are as in Table1.
a Total LAI (Larch LAI)

Table 2 .
Parameters used in the BIOME-BGC model.

Table 3 .
Climate biases applied into the sensitivity study.The bias for each season was determined as three times the standard deviation of the natural variability from the long-term climate data shown in Fig.2.

Table 4 .
Root mean square error (RMSE) between observed and simulated carbon and water fluxes at a daily time scale.

Table 5 .
The sensitivity of annual GPP, RE, and NEE to changes in weather patterns of air temperature, precipitation, VPD, and solar radiation.The numbers in the table are the relative percentage as compared to the control experiment.The applied climate bias in each season was determined as 3 standard deviations of the natural variability (Table3).The underscores in the table represent the values that show a high sensitivity for changing GPP and RE (more than 10%) and NEE (more than 20%) compared to the control experiment.