Modeling the impact of drought on canopy carbon and water fluxes for a subtropical evergreen coniferous plantation in southern China through parameter optimization using an ensemble Kalman filter

Soil and atmospheric water deficits have significant influences on CO 2 and energy exchanges between the atmosphere and terrestrial ecosystems. Model parameterization significantly affects the ability of a model to simulate carbon, water, and energy fluxes. In this study, an ensemble Kalman filter (EnKF) and observations of gross primary productivity (GPP) and latent heat (LE) fluxes were used to optimize model parameters significantly affecting the calculation of these fluxes for a subtropical coniferous plantation in southeastern China. The optimized parameters include the maximum carboxylation rate ( Vcmax), the slope in the modified Ball-Berry model ( M) and the coefficient determining the sensitivity of stomatal conductance to atmospheric water vapor deficit ( D0). OptimizedVcmax andM showed larger variations thanD0. Seasonal variations of Vcmax andMwere more pronounced than the variations between the two years. Vcmax andM were associated with soil water content (SWC). During dry periods, SWC at the 20 cm depth explained 61% and 64% of variations of Vcmax andM, respectively. EnKF parameter optimization improved the simulations of GPP, LE and SH, mainly during dry periods. After parameter optimization using EnKF, the variations of GPP, LE and SH explained by the model increased by 1% to 4% at half-hourly steps and by 3% to 5% at daily time steps. Further efforts are needed to differentiate the real causes of parameter variations and improve the ability of models to describe the change of stomatal conductance with net photosynthesis rate and the sensitivity of photosynthesis capacity to soil water stress under different environmental conditions. Correspondence to: W. Ju (juweimin@nju.edu.cn)


Introduction
Carbon (C) and water fluxes between the atmosphere and terrestrial ecosystems are interactively linked at various temporal and spatial scales (Rodrigues-Iturbe, 2000).Two major components of the terrestrial C cycle (C fixation through photosynthesis and C release through heterotrophic respiration) are strongly regulated by soil water content (SWC) (Barr et al., 2006).Consequently, SWC is a key determinant of net ecosystem productivity (NEP) of a terrestrial ecosystem (Grainer et al, 2003).Photosynthesis peaks when SWC is near the field capacity and decreases with the departure of SWC from this optimum.When soil or atmospheric water is insufficient, leaf stomatal conductance reduces and correspondingly CO 2 fixation is limited (Dang et al., 1997;Hogg et al., 2000;Grant et al., 2006a).Extensive studies have demonstrated the influence of SWC on C fixation at leaf (Infante et al., 1999;Lawlor, 1995) and canopy (Reichstein et al., 2002;Rambal et al., 2003) levels.For three Mediterranean evergreen forests, drought caused lightsaturated ecosystem gross C fixation and day-time averaged canopy conductance to decline by up to 90% and ecosystem water-use efficiency relative to the gross C fixation to decrease (Reichstein et al., 2002).The heat and drought hitting Europe in 2003 caused a 30% reduction in the continental gross primary productivity (GPP) (Ciais et al., 2005).The sensitivity of GPP to drought depends on vegetation species.In the boreal region of Canada, the decreasing degree of GPP caused by drought during 2001 to 2003 was larger for Aspen than for Jack Pine and Black Spruce (Kljun et al., 2005).The heat and drought in 2003 in Europe caused larger reductions in GPP in temperate deciduous beech and northern Mediterranean forests than in other ecosystems (Ciais et al., 2005;Granier et al., 2007).
Published by Copernicus Publications on behalf of the European Geosciences Union.
W. Ju et al.: Impact of drought on canopy carbon and water fluxes Models are useful tools to estimate regional/global terrestrial C budget.However, large discrepancy exists in outputs by models due to differences in model structure, input data, and parameters.The reliability of a model to simulate terrestrial C and energy fluxes depends, to a great extent, on its ability to simulate SWC dynamics and the influence of SWC on C sequestration.Performance of models in simulating the decline of C sequestration caused by soil and atmospheric dryness is quite different (Grant et al., 2006b).Accurate estimation of model parameters is critically important for modeling surface fluxes (Wolf et al., 2006).However, not all model parameters are measurable, especially at canopy level.And some parameters are not temporally constant.Properly determining parameter values is always a big challenge for modelers.
The availability of long-term measurements of CO 2 , water and energy fluxes using the eddy covariance (EC) technique in different ecosystems provides opportunities to optimize model parameters to which simulated CO 2 , water, and energy fluxes are very sensitive.Various parameter optimization methods have been successfully adopted in combination with EC measurements to optimize model parameters and to analyze the seasonality of parameters.Parameter optimization can be implemented by minimizing a cost function constructed on the basis of one or several criteria quantifying the differences between simulated fluxes and corresponding measurements (Aakto et al., 2004;Owen et al., 2007;Reichstein et al., 2003;Wang et al., 2007;Wolf et al., 2006) or by using the ensemble Kalman filer based approach (Mo et al., 2008).The cost function can be resolved at different temporal intervals to identify the seasonal variations of parameters.Through inverting a mechanistic model every five days, Reichstein et al. (2003) found that both active leaf area and photosynthetic capacity (V c max ) declined during the drought period for three Mediterranean forest ecosystems.With the inclusion of drought-dependent parameter variations, the model's ability to simulate hourly gross CO 2 uptake and transpiration was greatly improved.In a pristine grassland-forb steppe, inverted V c max and m (the Ball-Berry coefficient) also had large seasonal variations (Wolf et al., 2006).Recently, Mo et al. (2008) studied the trajectories of model parameters for an old aspen stand in Canada at daily time steps using an ensemble Kalman filter-based algorithm with EC measurements of GPP, ecosystem respiration and latent heat.Inversed parameter values (clumping index, V c max and m) showed distinct interannual and seasonal variations.Keenan et al. (2008) declared that for Mediterranean forest, changes in photosynthetic capacity not related to stomatal closure under water stressed conditions are very important for simulating surface carbon and water fluxes.
Water availability is an important ecological factor affecting C sequestration of terrestrial ecosystems in China (Yu et al., 2006).Long-term CO 2 flux measurements at 8 stations in China show a positive linear relationship between annual NEP and precipitation across biomes with different water availability conditions.In tropical and subtropical regions of China with monsoon climates, seasonal and interannual variability of precipitation are considerable.In these regions, although annual precipitation is normally plentiful, seasonal droughts occur frequently and influence terrestrial C sequestration.The decease in terrestrial C sequestration under seasonal drought conditions is determined by the accumulation of soil water deficits and a co-occurrence of high temperatures (Sun et al., 2006).On the basis of modeling analysis, Gu et al. (2008) found that drought caused NEP of a subtropical coniferous plantation in southern China to decrease by 63% and 47% in 2003 and 2004, respectively.At this site, C fixation was more sensitive to drought than ecosystem respiration.Soil drought has a much larger influence on NEP than atmospheric drought.Therefore, it is of importance for reliable estimation of terrestrial C balance in China to study the impact on C sequestration by soil and atmospheric water deficits.Although some observational and modeling research have been conducted to study the response of C sequestration by forests to drought in subtropical region of China (Sun et al., 2006;Gu et al., 2008), drought driven changes of model parameters have not been studied yet.
The objectives of this study are: (1) to study the interannual and seasonal variations of model parameters significantly affecting the calculations of canopy C and water exchange; (2) to analyze the linkage between variations of parameters and SWC; (3) to investigate the improvement in simulated canopy C and water fluxes by the inclusion of drought-dependent parameters.To achieve these objectives, the BEPS (Boreal Ecosystem Productivity Simulator) model (Ju et al., 2006) was optimized using eddy covariance measurements of GPP and latent heat fluxes (LE) taken over a subtropical coniferous plantation at QianYanZhou ecological station (QYZ), located in southeastern China.

Site description
The QYZ ecological station (115 • 01 13 E, 26 • 44 48 N) is located at the southwest portion of Jiangxi province (Fig. 1).This site belongs to the typical red earth hilly region in the mid-subtropical monsoon landscape zone of southeastern China.The mean annual temperature is 17.9 • C and annual precipitation is 1542.4mm during 1985-2002.The vegetation is mainly evergreen coniferous trees planted after 1985, including Pinus massoniana, Pinus elliottii engelm, and Cunninghama Lanceolata.An above-canopy flux system mounted in a 39.6 m tower has been in operation to record half-hourly exchanges of CO 2 , H 2 O and energy fluxes between the ecosystem and atmosphere and meteorological conditions since September, 2002 (Liu et al., 2006).This site is a typical plantation ecosystem in Jiangxi province, in which vegetation restoration (mainly reforestation and afforestation) has been widely and successfully practiced in the past three decades.Therefore, parameters optimized at this site can be used for investigating the effect of vegetation restoration on C sequestration in this region.

Data used
Half-hourly measurements of air temperature, incoming solar radiation, relative humidity, wind speed and precipitation during 2003 and 2004 were used for driving the BEPS model.The forest at QYZ is an evergreen coniferous plantation.The seasonal variability of leaf area index (LAI) is relative small (Huang et al., 2007).LAI was estimated from photosynthetically active radiation (PAR) measured above and below the canopy in clear days (Baldocchi et al., 1994).These estimated LAI values were averaged to produce monthly mean LAI, which was used to drive the model.Figure 2 shows the temporal trend of LAI, varying from 3.2 to 3.9 in two years.The maximum of LAI is larger in 2004 than in 2003.
Measurements of latent heat (LE) and CO 2 fluxes using the EC method (Liu et al., 2006) are used for parameter optimization.The missed NEP data were filled using the Michaelis-Menten equation (Liu et al., 2006).Missed LE data were filled using the Penman-Monteith equation.In order to reduce the influence of observation noises on parameter optimization, the flux data were checked in two steps.The methodology of Song et al. (2006) was first used to exclude abnormal records of half-hourly NEP and LE data.Then, if the summation of LE and sensible heat fluxes are below 70% or above 125% of the measured above-canopy net energy flux, the LE data was marked as abnormal values.GPP was derived from measured NEP, air temperature, and SWC (Wen et al., 2006).The averages of measured GPP and LE for the period between 9:00 and 17:00 were calculated for parameter optimization.

Model description
The BEPS model used in this study includes a module for photosynthesis calculation (Chen et al., 1999), a soil biogeochemical module for soil C, N and heterotrophic respiration simulations following CENTURY (Parton et al., 1993), and a land surface scheme for the computations of energy balance, sensible and latent heat fluxes, soil water and soil temperature status (Ju et al., 2006).GPP is calculated by scaling leaflevel Farquhar's biochemical model (Farquhar et al., 1980) up to canopy-level implemented through a sunlit and shaded leaves separation algorithm (Chen et al., 1999).Fluxes of C, water, and energy are separately simulated for sunlit and shaded leaves.Details about the model are described in Ju et al. (2006).We only highlight processes directly related to parameter optimization here.
GPP is calculated as: whereA net,sunlit and A net,shaded are net photosynthesis rates of sunlit and shaded leaves (mol m −2 s −1 ), respectively, calculated using the leaf level Farquhar model (Farquhar et al., 1980;Chen et al., 1997); LAI sunlit and LAI shaded are the leaf area index of sunlit and shaded leaves, respectively.The bulk stomatal conductance of the sunlit and shaded leaves for water vapor (G s,i , mol m −2 s −1 ) is calculated using a modified version of Ball-Woodrow-Berry (Ball et al., 1987) empirical model following Leuning (1995): where G o,i is the residual conductance (mol m −2 s −1 ); f w is a parameter describing the sensitivity of G o,i to soil water availability, m is the Ball-Berry coefficient; C s,i is the CO 2 concentration at the leaf surface (µ mol mol −1 ); D s,i is the vapor pressure deficit at the leaf surface (kPa); and D 0 is an empirical parameter determining the sensitivity of stomatal conductance to vapor pressure deficit (kPa), and is the compensation point (µmol mol −1 ).
In Eq. ( 2), parameters f w and m work together to link stomatal conductance with the net photosynthesis rate.f w can be empirically determined by soil water content (Wang and Leuning, 1998), soil water potential (Zierl , 2003) or leaf water potential (Tuzet et al., 2003).In the optimization, f w and m are grouped together and Eq. ( 2) was rewritten as: where M is a parameter representing the sensitivity of stomatal conductance to net photosynthesis rate, named as the slope of the modified Ball-Berry model in the flowing sections.LE is simulated as: where T r is the transpiration rate from canopy (kg m −2 s −1 ); E i and E g are evaporation rates of intercepted water from canopy and ground surface (kg m −2 s −1 ), respectively; and λ is the latent heat of water vaporization.Sublimation is not considered at this subtropical forest site.Canopy transpiration is calculated as the summation of transpired loss from sunlit and shaded leaves.Following Wang and Leuning (1998), transpiration from sunlit (i = 0) and shaded (i = 1) leaves is computed as: where ρ is the air density (kg m −3 ); C p is the specific heat of air at constant temperature (=1010 J kg −1 • C −1 ); γ the psychrometric constant (kPa • C −1 ); D a is the atmospheric vapor pressure deficit (kPa); is the rate of change of the saturated vapor pressure with temperature (kPa • C −1 ); T s,i and T a are temperatures at leaf surface and of air ( • C), respectively; r a and r b are boundary layer and aerodynamic resistance (s m −1 ), respectively.

Parameter optimization
Model parameter optimization was implemented using a dual state-parameter estimation method based on the ensemble Kalman filter (Moradkani et al., 2005;Mo et al., 2008).The parameter ensembles are updated as: where θ i+ t+1 is the updated parameter ensemble member at time step t+1; θ i− t+1 is the parameter sample generated with an assumption that model parameters follow a random walk in the space prescribed (Moradkhani et al., 2005); K θ t+1 is the Kalman gain matrix for correcting the parameter trajectories; and y i t+1 and ŷi t+1 are the ensemble members of observations and predictions, respectively Parameter sample θ i− t+1 is made as: In above Eq.( 6), the Kalman gain matrix K θ t+1 is calculated as: where θy t+1 is the cross covariance between parameter ensemble and prediction ensemble; yy t+1 is the error covariance matrix of the prediction ŷi t+1 ; and y t+1 is the error covariance of observation.

Parameters optimized and ensemble
Through sensitivity analysis, three parameters (i.e., V cmax , M, and D 0 ) are selected for optimization.These parameters have significant impacts on simulated GPP and LE.Although the electron transport rate (J max ) is also an important parameter for calculating GPP and LE, it was calculated according to V cmax in the BEPS model (Chen et al., 1999).Therefore, J max was not optimized in this study.The ranges and standard deviations of these parameters are listed in Table 1.The ensemble members of predicted GPP and LE are generated by disturbing these parameters using Eq. ( 7) with an assumption that data used to force the model are perfect.
Estimation of observational errors is important for EnKF data assimilation since the observational errors determine the degree to which the simulated results are to be corrected (Mo et al., 2008).Studies show that the observation errors of flux data are related to flux magnitudes (Lasslop et al., 2008, Richardson et al., 2008).According to the previous study on the uncertainties of observed flux data conducted by Liu et al. (2009) for this site, the uncertainties of GPP and LE are set as 15% of their daily averages and assumed to be independent of each other.The ensemble size is set to 200 to ensure the correct estimate of cross covariance between parameter and prediction ensembles and that of the prediction error covariance.Parameter optimization is conducted in days when both GPP and LE measurements are available.The parameters are updated by adding up corrections calculated using Eq. ( 6).Then the new parameters are used to simulate fluxes at half-hourly time steps in the following days until next parameter optimization is conducted.

Interannual and seasonal variations of the parameters
The seasonal variations of three optimized parameters were analyzed with 10-day averages since there are certain daily fluctuations of retrieved parameters.Figure 3 shows the seasonal variations of Vc max in years 2003 and 2004.This parameter was in the range from 14.4 to 29.3 µmol m −2 s −1 in 2003 and from 13.3 to 28.6 µmol m −2 s −1 in 2004, respectively.
The averages of Vc max were 19.2 and 20.1 µmol m −2 s −1 in 2003 and 2004, respectively.In 2003, this parameter gradually increased to a peak value of 29.3 µmol m −2 s −1 in the middle of June and then quickly decreased due to the impact of severe summer drought.It started to recover with the increase of precipitation in mid August and approached the second peak of 23.9 µmol m −2 s −1 in the middle of September.Then this parameter declined again and showed a relatively small variation after middle October.The seasonality of V c max in 2004 was very similar to that in 2003.The summer decrease of V c max was smaller in 2004 than in 2003.Starting from late September, this parameter was larger in 2004 than in 2003.
The parameter of V c max represents the foliage photosynthetic capacity at a reference temperature (25 • C in this study) and is believed to be linearly related with leaf nitrogen content (Chen et al., 1999).It is also affected by other factors including photosynthetic enzyme capacity, leaf mass area ratio, and leaf age, soil water stress, acclimation to seasonal change of environmental conditions (Wilson et al., 2006;Mo et al., 2008;Keenan et al., 2009).Many studies have demonstrated the seasonal variations of V c max (Reichstein et al., 2003;Wolf et al., 2006;Wang et al., 2007;Mo et al., 2008).However, the reported seasonal patterns of V c max varied considerably in different ecosystems.The seasonal changes of V c max identified in this study are similar to those in three Mediterranean forests reported by Reichstein et al (2003).In both cases, V c max showed increasing trends in the spring and early summer and a noticeable decline during a summer drought period.
The slope in the Ball-Berry model (M) ranged from 7.5 to 13.4 in 2003 and from 7.0 to 12.9 in 2004, respectively (Fig. 4).It averaged 10.3 and 10.7 in these two years.In 2003, it showed an increasing trend prior to the end of March and a decreasing trend after April.It dropped during the summer drought and recovered when the soil was getting wet after the middle of August.This parameter reached the second peak of 11.7 in mid September.Then, it continuously decreased as the soil was getting dry again.The seasonal pattern of M in 2004 is different from that in 2003.In 2004, this parameter gradually increased from mid March to early September and then slowly decreased.It did not decrease during the summer of 2004.
The seasonal and interannual variations of parameter M have been confirmed in many studies (Sala and Tenhunen, 1996;Lai et al., 2000;Wolf et al., 2006;Mo et al., 2008).There is still no consistent conclusion on the seasonal pattern of this parameter and its association with SWC.Sala and Tenhunen (1996) used single values of LAI and V c max for the entire growing season to invert M and found that this parameter changed dramatically.Lai et al. (2000) found that M was associated with SWC.Wolf et al. (2006) demonstrated that M could vary seasonally in the range from 15.4 to 23.8 and the association between M and SWC was not strong.Mo et al. (2008) indicated that M was smaller in dry years than in wet years for an aspen forest in the southern study area of BOREAS, Canada.The variations of M between two years in this study were similar to the results reported by Mo et al. (2008).Both of studies found that M was smaller in dry than in wet years.
The annual averages of D 0 were almost identical in 2003 and 2004, equal to 1.47.This value is close to those used in many modeling studies (Wang and Leuning, 1998;Arrain et al., 2003;Ju et al., 2006).This parameter showed similar seasonal patterns in 2003 and 2004.The seasonal variations of this parameter were relatively small, ranging from 1.35 to 1.64 kPa in 2003 and from 1.33 to 1.62 kPa in 2004 (Fig. 5), respectively.During the summer drought period in 2003, it decreased slightly and had smaller values than those in the same period of 2004.D 0 was slightly below the average during the periods from April to early August in 2003 and from mid February to late May in 2004.

The relationship between the parameters and soil water content
As expected, parameters V c max and M showed a certain degree of association with SWC (Fig. 6).These two parameters started to decline when SWC decreased.They are more tightly linked with SWC in deep layers than with SWC in the surface layer.V c max kept increasing when SWC at 5 cm depth started to decrease and SWC at depths of 20 cm and 40 cm remained above the field capacity.This confirms www.biogeosciences.net/7/845/2010/Biogeosciences, 7, 845-857, 2010  the hypothesis that forest is able to uptake deep water to maintain normal photosynthesis and transpiration rates during the initial stage of drought.With the continuation of drought, deep soil water is consumed and the stomatal conductance decreases (Ju et al., 2006).Consequently, the photosynthetic capacity (parameterized as V c max in the model) showed a considerable reduction (Reichstei et al., 2003).
V c max showed a stronger response to decrease in SWC in mid summer of 2003 than in mid summer of 2004 due to cooccurrence of drought and high temperature in 2003 (Sun et al., 2006).V c max and SWC at depths of 20 cm and 40 cm almost started to decrease simultaneously.However, the recovering rate of V c max was slower than that of SWC, indicating drought has a lagged effect on the photosynthetic capacity.In years 2003 and 2004, the lags with maximum correlation with V c max were about 20 days and 10 days for SWC at depths of 20 cm and 40 cm, respectively.Parameter M decreased during the summer drought period in 2003.This decrease of M did not occur during the summer of 2004 although V c max showed a moderate decrease because SWC of deep layers were relatively high.
When SWC at 20 cm depth was below the field capacity, both M and V c max linearly increased with SWC.If SWC was above the field capacity, its change had no detectable influence on M and V c max (Fig. 7).This is similar to the findings of the impact of SWC on M in previous studies (Lai et al., 2000;Wolf et al., 2006;Kennan et al., 2009) and on V c max by Reichstein et al. (2003) and Wang et al. (2007).SWC alone explained 61% and 64% of the variations in V c max and M when SWC was below the field capacity.V c max and M showed very large variations when SWC was above the field capacity, indicating that SWC was the key determinant of V c max and M when SWC was low, causing stress on forest growth.However, when SWC was high, other factors may also have significant influences on them.

Performance of the optimized BEPS model in flux simulation
Because fluxes of GPP, LE and sensible heat flux (SH) were predicted with the corrected parameters on the previous day, the efficiency of model parameter optimization in improving flux estimation was evaluated by comparing measured and simulated fluxes.Figure 8 shows the comparisons for half-hourly GPP, LE and SH.The simulated fluxes are the averages of 200 ensemble predictions.With the optimized parameters, the BEPS model explained 83%, 81% and 83% of variations of measured half-hourly LE, SH and GPP in 2003, respectively.In 2004, the optimized model captured 83%, 79% and 83% of variations in measurements of LE, SH and GPP.The model tends to overestimate smaller values of LE and GPP and underestimate their large values.In contrast, the small (or large) values of SH were underestimated (or overestimated) by the model, respectively.These results are significant improvements over the results simulated with stomatal conductance calculated using Eq. ( 2) with calibrated parameter set (D 0 =1.58 kPa, m = 11, V c max =25 µmol m −2 s −1 ) (Table 2).The biggest improvement was achieved for GPP, followed by LE and SH.The RMSE (root mean square error) values of GPP with EnKF parameter optimization decreased to 0.17 and 0.19 µmol m −2 s −1 in 2003and 2004, respectively, from 0.21 Biogeosciences, 7, 845-857, 2010 www.biogeosciences.net/7/845/2010/At daily steps, the optimized model was able to explain 90% and 91% variations of daily GPP in 2003 and 2004, respectively (Fig. 9).These percentages for the simulations using the calibrated parameter set were 86% and 89%.The improvement for daily GPP simulation was mainly in the periods of the summer and autumn droughts of 2003 and 2004.The improvement in simulated GPP after EnKF parameter optimization was smaller in 2004 than in 2003.With the optimized parameters, the model greatly reduced the overestimation of GPP during the drought periods (Fig. 9).
The improvement in the simulation by the EnKF parameter optimization was greater for LE than GPP (Figs. 9 and 10).93% and 90% of variations of day time LE were captured by the BEPS model optimized using EnKF in 2003 and 2004, respectively while corresponding values were 89% and 85% for simulations using the calibrated parameter set in 2003 and 2004.The improvements in LE simulations were also mainly during the drought periods in 2003 and 2004.During the serious drought period from day 171 to day 232 in 2003, the RMSE of simulated LE was 31 W m −2 using the optimized parameters while it was 43 W m −2 using the calibrated parameters.During the drought period from day 188 to day 221 in 2004, the RMSE of simulated LE was 41 W m −2 using the optimized parameters while it was 56 W m −2 using the calibrated parameters.
The model using the optimized parameters was unable to capture some peaks of GPP and LE in the midday since the adjustments of parameters was carried out using daily values of observation GPP and LE.

Discussions
Parameters in the ecosystem models are temporally and spatially interlinked in nonlinear forms.Parameter values derived from measurements taken at leaf or canopy scales may lead to incorrect predictions at other scales (Mo et al., 2008).Eddy covariance measurements represent the information at 100-1000 m footprint distance, similar to the spatial resolution of many hydrological and ecological models.Therefore, parameter values optimized using eddy covariance measurement will improve the simulations at regional scales.www.biogeosciences.net/7/845/2010/Biogeosciences, 7, 845-857, 2010  Many studies have demonstrated the applicability of eddy covariance measurements to optimizing model parameters (Reichstein et al., 2003;Wolf et al., 2006;Wang et al., 2007;Santaren et al., 2007;Mo et al., 2008).EnKF sequential data assimilation is an effective technique for model parameter optimization, which is able to track the seasonal variations of the parameters representing the photophynthetic capacity indicated by V c max and the sensitivity of stomotal conductance to net photosynthesis rate indicated by M. The estimation of observation errors is critical to the EnKF parameter optimization algorithm.If larger observation errors are used in the calculation of Kalman gain matrix, the trajectories of model parameters will be corrected to a smaller degree.Using the daily-differencing approach, Liu et al. (2009) found that relative errors of observed GPP decreased with temporal scales and the overall uncertainties of observed GPP at half-hourly and daily time scales are about 10% during 2003 to 2005 at this site.Richardson et al.(2008) and Lasslop et al.(2008) proposed practical methods to estimate observational errors of flux data and pointed out that the observation errors of flux data increased with flux magnitudes.In this study, the errors of GPP and LE were set as 15% of their daily values, which are in the range of reported uncertainties of GPP and LE measurements (Santaren et al., 2007;Richardson et al., 2008) and close to the numbers reported by Liu et al. (2009).The assumption ignores the effects of meteorological conditions on the proportion of observational error of flux data to flux magnitudes, possibly inducing some uncertainties in optimized parameters.The application of estimated observational errors of flux data for each day will increase the confidence of optimized parameters.The construction of the second eddy covariance tower in July, 2008 at this site will make it possible to estimate uncertainties of GPP and LE measurements more reasonably.
Current conclusions on the statistical properties of random flux errors are not inconsistent.Hollinger and Richardson (2005) indicated that flux measurement errors are more close to a double exponential distribution.Liu et al. (2009) declared a similar conclusion about the statistical distribution of random errors of observed fluxes for the QYZ site.Recently, Lasslop et al. (2008) indicated that the double exponential distribution of random errors of observed errors is mostly due to the superposition of normal distributions with varying standard deviation.In this current study, the observation error is assumed to be Gaussian distributed and independent between GPP and LE.The effects of this assumption about the distribution of observation errors on optimized parameters need further investigation.
Compared with the global optimization algorithms, which use eddy covariance measurements in long time periods to invert parameters, the EnKF optimization method updates the model parameters using measurements in one day.Although   this strategy is able to trace day to day variations of parameters, it makes the optimization results sensitive to the quality of observation.One or several abnormal observational values in one day might cause the parameters to vary unrealistically.Therefore, careful quality control of observation data is of importance prior to the implementation of the EnKF parameter optimization scheme.Since the daily averages of observed and simulated GPP and LE are used to correct the parameters, the simulation accuracy is more improved at daily time steps than at half-hourly time steps.This is because in the calculations of daily values, the underestimation/overestimation of large GPP or LE values will be possibly compensated by overestimation/underestimation of small GPP or LE values.The model with optimized parameters still tends to underestimate large GPP and LE fluxes at noon.Similar to findings in other studies, considerable seasonal variations of V c max and M are derived through assimilating eddy covariance.They also show noticeable changes between the two years used in our study.There are many possible mechanisms causing the changes of these parameters.In addition to water stress, changes in photosynthetic enzyme capacities, foliage N concentration status, leaf mass per area, leaf age and the acclimation of forest to environmental conditions might cause the seasonal variations of V c max (Wilson et al., 2000;Reichstein et al., 2003;Mo et al., 2008).Carlhais et al. ( 2008) indicated that model structure and estimation method can also cause the temporal variations of optimized parameters.In this study, the model structure and optimization method are assumed to be perfect.The robustness of this assumption needs further study.
Development of algorithms to quantify the variations of these parameters (V c max and M) is of importance for improving the simulations of carbon and water fluxes at various scales since they are used in models to represent the photosynthetic capacity and the sensitivity of stomatal conductance to net photosynthetic rates.Previous studies indicated that the allowance of the change in photosynthetic capacity (V c max ) with soil water stress was required to improve the agreement between simulated and observed fluxes during dry periods (Reichstein et al., 2003;Wang et al., 2007).Current study found that both responses of the sensitivity of stomatal conductance to net photosyntheis rates (indicated by the reduction of M) and photosynthetic capacity to soil water stress (indicated by the reduction of V c max ) are necessary for reliable simulations of GPP and LE for this subtropical coniferous forest ecosystem.This is similar to the findings for the Mediterranean forest ecosystem reported by Keenan et al. (2009).
It was found that SWC is a good predictor for V c max and M during dry periods for this coniferous plantation ecosystem in subtropical zone of southeastern China.During wet periods, V c max and M showed large variations around the averages (Fig. 7).The simultaneous observations of nitrogen content, V c max , M, photosynthesis rate are required to understand the mechanisms underlying the variations of V cmax and M and to develop methodologies to describe the variations of V cmax and M for bettering the simulations of surface fluxes.

Conclusions
From this study, following conclusions can be drawn: 1.The EnKF data assimilation is an effective technique for model parameter optimization, which is able to track the seasonal and interannual variations of the parameters.
The optimized parameter values are sensitive to the estimation of observation errors and quality of observation data.
2. Optimized V c max and M showed considerable seasonal variations and also changed noticeably between two years.The seasonal variations of parameters were larger than their variations between the two years.The variations of V c max and M indicate that soil water stress reduced the sensitivity of stomotal conductance to net photosynthesis rates and photosynthesis capacity.SWC is a good predictor of V c max and M during dry periods with SWC of 20 cm depth below the field capacity.SWC of 20 cm depth alone explained 61% and 64% variations of V c max and M, respectively.During wet periods, the dependence of these two parameters on SWC were not significant and other factors (e.g., photosynthetic enzyme capacities, foliage N concentration status, and leaf mass per area, acclimation of forest to environmental conditions and leaf aging) can cause these parameters to vary largely.
3. The inclusions of responses of both the sensitivity of stomatal conductance to net photosynthesis rate (M) and photosynthetic capacity (V c max ) to soil water stress improved the simulation of GPP, LE and SH under different soil water conditions.After parameter optimization using EnKF, the variations of GPP, LE and SH explained by the model increased by 1% to 4% at halfhourly steps and increased by 3% to 5% at daily time steps.
4. In this study, the effects of model structure and optimization method (causing stochastic variations) on the optimization of parameters were not investigated.In addition, the observational errors can be split into systematic and random components.Only the random components of observational errors were currently included.Further efforts are required to consider these effects of model structure and optimization methods on parameter optimization to identify the real mechanisms underlying the temporal variations of parameters.

Fig. 1 .
Fig. 1.The location of the study site.Green colored areas are forests.The red solid circle represent the QianYanZhou Ecological station.

Fig. 2 .
Fig. 2. Monthly LAI derived from photosynthetically active radiation measured above and below the canopy in years 2003 and 2004.Error bars represent standard deviations of estimated LAI in each month.

Fig. 4 .
Fig. 4. Seasonal variations of optimized M in 2003 (left) and 2004 (right).Error bars represent standard deviations of optimized M within 10 days.

Fig. 6 .
Fig. 6.Seasonal variations of optimized parameters of V c max and M and measured soil water content at depths of 5 cm, 20 cm and 40 cm.The left panel is for year 2003 and the right one is for year 2004.

Fig. 7 .
Fig. 7. Association of optimized Vc max and M with measured relative soil water content (fraction of field capacity) at 20 cm depth in 2003 and 2004.Values shown are 10-day averages.Curves are linear functions with upper threshold, fitted to unfilled symbols.Also indicated are the R 2 values of fitted linear functions.Error bars are the standard deviations of optimized parameters within 10 days.

Fig. 8 .
Fig. 8.Comparison between observed and simulated half-hourly fluxes of GPP, LE, and SH.The grey dashed line is the 1:1 line and dark solid one is the regression line.The left panel is for 2003 and right one is for 2004.

Fig. 9 .
Fig. 9. Residuals of simulated daily GPP using parameters optimized by EnKF and parameters determined through model calibration for years 2003 (left) and 2004 (right).

Fig. 10 .
Fig. 10.Residuals of simulated daily mean LE using parameters optimized by EnKF and parameters determined through model calibration for years 2003 (left) and 2004 (right).

Table 1 .
The standard deviations and ranges of three parameters optimized