Variability of phenology and fluxes of water and carbon with observed and simulated soil moisture in the Ent Terrestrial Biosphere Model ( Ent TBM version 1 . 0 . 1 . 0 . 0 )

The Ent Terrestrial Biosphere Model (Ent TBM) is a mixed-canopy dynamic global vegetation model developed specifically for coupling with land surface hydrology and general circulation models (GCMs). This study describes the leaf phenology submodel implemented in the Ent TBM version 1.0.1.0.0 coupled to the carbon allocation scheme of the Ecosystem Demography (ED) model. The phenology submodel adopts a combination of responses to temperature (growing degree days and frost hardening), soil moisture (linearity of stress with relative saturation) and radiation (light length). Growth of leaves, sapwood, fine roots, stem wood and coarse roots is updated on a daily basis. We evaluate the performance in reproducing observed leaf seasonal growth as well as water and carbon fluxes for four plant functional types at five Fluxnet sites, with both observed and prognostic hydrology, and observed and prognostic seasonal leaf area index. The phenology submodel is able to capture the timing and magnitude of leaf-out and senescence for temperate broadleaf deciduous forest (Harvard Forest and Morgan– Monroe State Forest, US), C3 annual grassland (Vaira Ranch, US) and California oak savanna (Tonzi Ranch, US). For evergreen needleleaf forest (Hyytiäla, Finland), the phenology submodel captures the effect of frost hardening of photosynthetic capacity on seasonal fluxes and leaf area. We address the importance of customizing parameter sets of vegetation soil moisture stress response to the particular land surface hydrology scheme. We identify model deficiencies that reveal important dynamics and parameter needs.


Introduction
Phenological timing remains a major weakness of land surface dynamic global vegetation models (DGVMs) that are coupled to general circulation models (GCMs) and a primary cause of uncertainty in predicting the trajectory of global atmospheric CO 2 (Friedlingstein et al., 2006(Friedlingstein et al., , 2014)).Seasonal variation of vegetation foliage, i.e., leaf phenology, determines the timing and duration of the photosynthetically active canopy, influencing stomatal activity, surface albedo and surface roughness (Jolly and Running, 2004).Thus, it plays a crucial role in the exchange of water, energy and carbon between land and the overlying atmosphere.Numerous observations show that the interannual variability of transpiration and gross primary productivity is associated with timings of leaf-out and leaf senescence across ecosystem types (Goulden et al., 1996).Light-controlled leaf phenology is suggested as a key controlling factor responsible for increasing carbon and water fluxes from land to the atmosphere during the dry season in the Amazon rainforests (Hutyra et al., 2007;Kim et al., 2012).Phenology is also tightly connected to other ecosystem processes, exerting strong controls on the amount of assimilated carbon that is subsequently utilized for plant growth and reproduction.Kramer et al. (2000) showed that phenology could have effects on the species composition of temperate-zone deciduous forests and the geographical distribution of species since difference in phenological response leads to difference in light availability and therefore growth in mixed species stands.Given the strong interactions Published by Copernicus Publications on behalf of the European Geosciences Union.

Y. Kim et al.: Variability of phenology and fluxes of water and carbon in Ent TBM
between phenology and other land surface and ecosystem processes, phenology affects both weather and climate.Seasonal variation in vegetation characteristics have been shown to significantly influence summer precipitation and temperature in the US (Dirmeyer, 1994;Xue et al., 1996) and enhance or weaken the feedbacks between soil moisture and precipitation in the continental interior of North America depending on soil moisture conditions and season (Kim and Wang, 2007).Levis and Bonan (2004) demonstrated that the coupling between phenology and the atmosphere is critical for models to capture seasonal weather evolution.Tightly linked to phenology, plant carbon allocation, that is, distribution of assimilated carbon among the plant parts, also responds to environmental and climate conditions (such as increases in air temperature, changes in precipitation patterns and elevated atmospheric CO 2 concentration).For example, Pumpanen et al. (2012) observed that root biomass and the rate of photosynthesis for silver birch, Norway spruce and Scots pine seedlings increase with higher soil temperature, yet a simultaneous increase in both photosynthesis and respiration rates results in no change in net CO 2 exchange and total seedling biomass.
Terrestrial biosphere models (TBMs) or dynamic global vegetation models have been developed and coupled to general circulation models (e.g., Foley et al., 1996;Cox, 2001;Sitch et al., 2003;Bonan and Levis, 2006;Dunne et al., 2013) to account for biophysical and biogeochemical processes and sometimes biogeography, allowing for the prediction of transient terrestrial ecosystem interactions with climate (Cramer et al., 2001;Friedlingstein et al., 2006).Thus, the active role of vegetation phenology can be incorporated into climate modeling.TBMs have been parameterized and evaluated on the basis of local, regional, or global scale studies.It has become common to evaluate the models at the individual field scale at sites with eddy flux measurements and detailed ground data (e.g., Delire and Foley, 1999;Arora and Boer, 2005;Krinner et al., 2005;Kucharik et al., 2006;Friend et al., 2007;Stöckli et al., 2008;Bonan et al., 2011).Still, parameterizations for vegetation processes (such as phenology and carbon allocation) implemented in TBMs are often limited to local-scale derivations due to the lack of high-quality global-scale observations of vegetation structure and function together with meteorological conditions and mechanistic understanding free of local effects.
Prognostic phenology models have been developed to predict the phenological response of vegetation to climate based on empirical evidence, as a mechanistic, process-based treatment is still not fully realizable with current understanding (Sala et al., 2012).The more commonly used climatic rulebased approach accounts for cues by temperature, soil moisture and day length cues to phenology, to predict leaf-on and leaf-off, with these controls often represented as cumulative functions of one or several climate variables that reach an empirically defined threshold (White et al., 1997).Another approach is based on plant carbon status (Bonan et al., 2003) and predicts leaf-out and senescence on the basis of potential positive carbon assimilation, which is in turn affected by temperature, moisture and sometimes nutrient conditions.
All of the above approaches require empirical parameterization of the responses to climate, and a model scheme that is independent of plant functional type (PFT) or geographical variation is still a research goal.Jolly et al. (2005) have proposed a very simple and promising bioclimatic growing season index (GSI) for phenology based on linear relations to minimum temperature, photoperiod and vapor pressure deficit (VPD, as a proxy for soil moisture), which seems to perform well compared to satellite observations at diverse sites.However, it performs less well for arid systems for which VPD may not be a good indicator of available deep soil moisture, and it is not able to capture any seasonal moisture or light sensitivity that has been observed in tropical evergreen forests (Stöckli et al., 2011).Forkel et al. (2014) adopted the concept of GSI but used the soil water availability instead of VPD for water limiting functions.Phenology depends not only on atmospheric water demand but also on water supply from soil moisture as Migliavacca et al. (2011) have shown that GSI performed better when using a soil moisture limiting function instead of the VPD limiting function.Recently, Caldararu et al. (2014) introduced a promising optimality approach based on the hypothesis that phenology is a strategy for optimal leaf area index, rather than explicit carbon exchange, driven by canopy-level demand for -and constrained by availability of -light and water, limited by leaf aging.They fitted the model to satellite observations of LAI (leaf area index) and demonstrated its capability to reproduce phenological patterns for different vegetation types over the globe within 8-16 days of observations.Top-down optimality approaches such as this may indeed be the smart way for global-scale models to capture the strategic behaviors inherent in phenology, in lieu of mechanistic understanding at the leaf or molecular level; the next step remains to couple them with explicit carbon exchange and allocation.
In this study, we perform a site-based model evaluation study for the Ent 1 Terrestrial Biosphere Model's (Ent TBM version 1.0.1.0.0 2 ) coupled phenology/growth schemes.This evaluation is a necessary task before introducing prognos-tic phenology into global simulations coupled with a GCM atmosphere in order to enable modeling of interactive phenology and climate.We do not offer yet a new paradigm, but the phenological timing schemes provide a synthesis of approaches in the literature to capture the full combination of climatological drivers thus far known to be essential for each type of phenology and introduce some new functional representations to do so.These are coupled to growth algorithms from the Ecosystem Demography (ED) model (Moorcroft et al., 2001) that account for both the geometric and mass allometry of PFTs.
In this paper, we describe the Ent TBM's phenology and allometry scheme coupled to the ED carbon allocation scheme and evaluate their performance at Fluxnet sites (Baldocchi et al., 2001), focusing on seasonal and interannual variations of LAI and carbon and water fluxes.We compare site simulations using both observed soil moisture and that modeled by a land surface hydrology model coupled to the Ent TBM.The phenology schemes synthesize several observational data sets, combining both climate responses and a carbon balance approach, described in detail below.Here we evaluate the performance for temperate broadleaf deciduous forest, C3 annual grassland, evergreen needleleaf forest and tree/grass savanna (mixed drought deciduous broadleaf and C3 annual grassland).Through these evaluations, we are interested in quantifying the accuracy of the current model at the site level, and we identify ecosystem processes needing further improvement, with regard to both plant growth dynamics and the representation of soil moisture.

Land surface model (LSM) of the NASA GISS GCM
The Ent TBM can be run with observed soil moisture and temperature and canopy temperature inferred from eddy flux measurements of sensible heat fluxes or, given precipitation and air temperature, it can obtain modeled soil moisture, temperature and canopy temperature if run coupled to a land surface hydrology model.For the coupled mode, we use the land model of the National Aeronautics and Space Administration (NASA) Goddard Institute for Space Studies (GISS) general circulation model (GCM) (Schmidt et al., 2006).The NASA GISS GCM land hydrology consists of six soil layers down to 3.5 m depth based on Rosenweig and Abramopoulos (1997), with updates described in Schmidt et al. (2006Schmidt et al. ( , 2014)).The LSM computes the fluxes of heat and water vapor to the atmosphere, and the energy balance of the soil and vegetation canopy.Surface runoff is calculated based on saturation and infiltration capacity of the upper soil layer.The underground runoff is computed according to a formulation of Abramopoulos et al. (1988), which takes into account the average slope and the density of underground sinks in the cell.
When running the Ent TBM coupled to the GISS LSM, soil physics parameters are taken from the land-surface-mapped data sets of the GISS LSM.

Ent Terrestrial Biosphere Model (Ent TBM)
The Ent TBM is a stand-alone model developed specifically for coupling the fluxes of water, energy, carbon and other trace gases between LSMs and GCMs.It is structured like the ED model (Moorcroft et al., 2001) for simulating competition in mixed canopies and disturbance dynamics by representing vertical canopy structure through ensemble cohorts of identical individuals and horizontal heterogeneity via subgrid patch communities.The specifications of canopy geometry and allometry of biomass pools are consistent with individual ellipsoidal crown geometry that is integrated with the coupled phenology/growth model.This paper presents simulations of seasonal variation in leaf area and mass and in fluxes of CO 2 , water vapor and sensible and latent heat of both transpiration and ground evaporation.
Figure 1 shows a conceptual diagram of the Ent model, and how it is coupled with a GCM (or offline meteorological forcings) and an LSM.Ent's biophysics modules operate at the physical time step of the GCM or LSM.The photosynthetic uptake of carbon utilizes the well-known photosynthesis model of Farquhar et al. (1980) and Farquhar and von Caemmerer (1982) coupled with the stomatal conductance model of Ball and Berry (Ball et al., 1987), while Ent uses its own cubic solution for these coupled equations.Canopy radiative transfer is optionally modeled as in Friend and Kiang (2006) for homogeneous canopies or as in Ni-Meister et al. (2010) and Yang et al. (2010) for clumped canopies.In this paper, in lieu of detailed site allometric and canopy structure data, we utilize the homogeneous canopy radiative transfer scheme.Carbon uptake is accumulated over a day so that carbon allocation to growth, phenological behavior and mortality are updated once per day.An individual plant has distinct biomass pools, including a "labile" or carbohydrate reserve pool into which photosynthetic uptake and re-translocated carbon are accumulated: "active" pools consisting of foliage, fine roots, a reproductive pool and, for woody plants, live sapwood; and "dead" pools consisting of dead stem wood and coarse roots.Autotrophic respiration is the sum of maintenance respiration as a function of biomass and temperature, "activity growth respiration" as a function of gross assimilation and tissue growth respiration as a function of amount of new growth.
Ent takes its meteorological drivers and hydrological balance at the grid-cell or catchment-zone scale of the LSM and subgrid heterogeneity is represented as dynamic patches of vegetation communities, comprised of cohorts of plants that are ensembles of identical individuals (patch and community dynamics are not part of this study).The biomass pools and geometry of an individual woody plant are illustrated in Fig. 2    to the grid-cell or catchment-zone level to couple with the atmosphere.Also, root density vertical profile distributions in Ent are used to calculate a depth-weighted average of soil moisture stress.These profiles are a modification of those in Rosenzweig and Abramopoulos (1997), with details given in Appendix A.
The Ent TBM is designed to support a flexible number of PFTs.A parameter set for 17 PFTs has been developed, as listed in Table 1; however, we note that only a subset of these PFTs is evaluated here according to data availability, and the others must be approximated from the available similar types and theoretical/empirical relations from the literature.Following the rationale first advocated by Defries et al. (1995) and adopted by all vegetation models since then to varying degrees, Ent's PFTs distinguish photosynthetic pathway (C3 and C4), phenological type (evergreen, cold deciduous and drought deciduous), leaf type (broadleaf and needleleaf), growth form (tree, shrub and herbaceous) and cultivation (herb crops).In addition, to better capture community dynamics in mixed canopies, if parameter sets are provided, Ent has the capability to distinguish early and late successional species through differences in leaf life span, following the approach of the ED model (Moorcroft et al., 2001), which is based on leaf physiological relations found in Reich et al. (2007).To capture total net carbon fluxes, the Ent TBM incorporates the code implementation of CASA' (Carnegie-Ames-Stanford Approach) from the Community Land Model 3.0 (CLM 3.0; Randerson et al., 2009;Doney et al., 2006; code kindly supplied by Jasmin John), which is based on Potter et al. (1993).For the Ent TBM, the CASA' temperature and soil moisture responses of respiration were replaced with functions derived from new fits to field data of Del Grosso et al. (2005).Details are provided in Appendix B.
As mentioned earlier, the Ent TBM can be run in several different modes of coupling: (1) a stand-alone mode when the meteorological (e.g., radiation, precipitation, air temperature, (2) a mode coupled with an LSM for prognostic soil moisture and temperature given meteorological forcings ("Ent-LSM"); and (3) a fully coupled mode with an atmospheric GCM.Ent-standalone and Ent-LSM modes can be used for site-specific simulations or regional/global simulations using observed meteorological and soil moisture data.
The Ent TBM can also be run with different levels of vegetation dynamics turned on or off.In a biophysics-only mode, canopy structure and leaf area are prescribed to simulate only fluxes of water vapor, carbon dioxide and other trace gases.In an "active biomass" phenology-only mode, canopy stem structure is prescribed and static, while seasonal leaf and fine root dynamics are prognostic, and carbon that would have been allocated to stem and coarse root growth instead is allocated to litter.In a phenology-woody growth mode, in addition to leaf phenology, stem and coarse root growth are also enabled.In an ecosystem-dynamics mode, mortality and disturbance ensure that plants cannot grow indefinitely and are subject to succession and cover change (ecosystem dynamics are not covered in this paper).

Plant growth submodel
The plant growth submodel integrates phenological timing and allocation of carbon to growth and litter fluxes (background litterfall and seasonal), and respiration fluxes are tied to tissue growth.The phenology scheme determines the phenological status of plants based on various environmental and climate rules studies, which determine budburst, frost hardening and senescence according to the phenological types of plants such as drought deciduousness and cold deciduous-ness.The carbon fixed over the course of each day from photosynthesis is accumulated and placed into a labile carbohydrate reserve pool.Carbon from the labile pool is then allocated once per day into different plant pools of foliage, sapwood, heartwood, fine root and coarse root as well as a reproductive pool according to empirical allometric relationships and leaf phenological status.In addition, tissue lost to background litter fluxes is replenished, and respiration fluxes are produced from growth of any tissue.A portion of litterfall is re-translocated back to the labile pool.
In the Ent TBM, the carbon allocation scheme takes a traditional approach of "static allocation", based on fixed allometric relationships between different pools, adopted from approaches of the ED models (Moorcroft et al., 2001;Medvigy et al., 2009).Appendix C provides the descriptions of the ED allocation scheme, which treats "active" and "dead" biomass pools as bulk sinks, with modifications for Ent.We identified some deficiencies of the ED allocation scheme and suggest future work for improvement in Sect. 5. Also note that Appendix D provides the biophysical, phenological and allocation parameter values used in this study.
Full prognostic growth entails growth of woody structure and the size of woody plants, which would require in addition full mortality and establishment dynamics so that there is no unlimited growth; these population and community dynamics will be presented in future papers.This study focuses on the "active biomass" performance of Ent given seasonal phenology, keeping woody structure static, allocating the amount that would have gone to growth to litterfall instead.

Phenology
The phenology scheme in the Ent TBM provides a synthesis and combines the climatic rule-based approach and carbon balance for deciduous plants to determine the timings and rates of leaf-out and leaf senescence by integrating several different modeling studies (Bonan et al., 2003;Botta et al., 2000;Foley et al., 1996;White et al., 1997).We present a diversity of PFTs, adding those with known behaviors that depart from common representations of cold, drought or light responses.While globally applicable parameterizations of climate rule-based phenology may still be elusive, where available in the literature, we draw from wide surveys that attempt to extrapolate to the global scale.
For deciduous plants, we use parameterizations by Botta et al. (2000).With growing degree day (GDD) and chilling requirement, they examined the possibility of extrapolating existing local models for leaf onset date to the global scale by retrieving leaf onset dates from the NOAA AVHRR (Advanced Very High Resolution Radiometer) satellite normalized difference vegetation index (NDVI).They identified appropriate leaf onset date models and estimated their parameters for each biome, which are implemented in other ecosystem models (Medvigy et al., 2008) Richardson et al. (2012), who conducted an intercomparison of phenology predictions of 11 TBMs (and three biophysics models with prescribed phenology) at five deciduous broadleaf and five evergreen needleleaf Fluxnet sites.They found that, for deciduous forests, the models consistently predicted an earlier onset of the growing season and later fall senescence than observed; meanwhile, most models underpredicted the magnitude of peak GDD sums, while those that explicitly or implicitly included a chilling requirement did relatively well in capturing the onset of LAI and GPP (gross primary production) for deciduous and evergreen forests, compared to simple temperature threshold schemes.For drought deciduous trees and grasses, we also make use of parameterizations of White et al. (1997) who developed a regional phenology model for the US, predicting timings of leaf onset and offset based on the satellite NDVI at the 20 km resolution.Their prediction errors are ∼ 1 week, and maximum expected errors are 10-14 days.
For evergreen vegetation, the Ent TBM includes frost hardening for boreal evergreen plants.The frost hardening (also called winter cold hardening) involves physiological changes to protect the plant from chilling injury and freezing injury, leading to a downgrading of leaf photosynthetic capacity as well as tissue turnover and respiration.Coniferous vegetation in the boreal zone has a clear annual cycle of photosynthetic activity, with photosynthesis low or zero during the winter, increasing during the spring, peaking during the summer and decreasing during the fall.While part of the cycle is due to direct responses to PAR (photosynthetically active radiation) and air temperature, the inherent photosynthetic capacity of needles also changes (Mäkelä et al., 2004).Therefore, the models that do not account for cold hardening and de-hardening will overpredict the uptake of carbon by photosynthesis for boreal systems during the late fall through early spring.This study implements a frost-hardening algorithm based on Hanninen and Kramer (2007), Mäkelä et al. (2006) and Repo et al. (1990), who developed a model of the frost hardiness of the stems of Scots pine seedlings.Below we describe the explicit model functions reflecting our choices based on the literature above.

Phenology model climate cue framework
In the Ent TBM, several "phenological factors", φ x , as well as physiological stress factors, β x , are calculated for seasonal environmental cues from various climate measures x.These include air and soil temperature history (cumulative number of growing degree days and of chilling days), day length and soil moisture.The phenological factors control the allocation of assimilated carbon, while the physiological stress factors affect the efficiency of carbon uptake, and all range from 0 to 1 on a daily basis.Different rules apply to the different PFTs, according to phenotype (woody plant cold deciduous, cd; drought deciduous, dd; evergreen, ev; tropical radiation phenology, tr; and cold deciduous herbs, c, whether annual or perennial).The phenological factor controls the timing and rate of carbon transfer between the labile and active carbon pools and hence the seasonal variation in LAI, fine roots and sapwood.
Furthermore, the Ent TBM determines "phenological status", Phenostatus p , where p is the phenotype, which identifies phenologically different seasons.For plants with seasonal leaf-out and senescence, Phenostatus p is 1 for the leafoff season, 2 for the leaf-up period, 3 for the peak foliage period and 4 for the senescent period.The trend in length of days (l day) is used to determine which season it is or, rather, which half of the year it is.If day length is decreasing, then it is the latter half of the year, and "fall" may be allowed to commence, depending on other climate variables of phenological factors.Below we itemize these variables and equations in the Ent phenology scheme.

Cold deciduous woody plants
During the winter, the phenological status of cold deciduous trees and shrubs, Phenostatus cd , is 1, for no foliage.Leafout (Phenostatus cd = 2) occurs once the cumulative number of GDDs exceeds its critical number (GDD crit ), which is determined with a function of cumulative number of chilling days (NCDs) (Botta et al., 2000).Following Kim and Wang (2005), the 10-day running average of air temperature (T 10 ) difference from the base temperature (T base ) of 5 • C is used to calculate GDD and NCD on a daily basis as follows: where t is time in days.GDD and NCD are reset to be zero at the beginning of the winter season (when Phenostatus cd switches from 4 to 1).The function for GDD crit is expressed as follows: where the constant values of GDD intercept , GDD slope and NCD multi are provided in Table 2.
Once leaf-out starts, trees take a number of degree days (GDD length ) to reach the phenologically unconstrained status (Foley et al., 1996).We introduce an approach to scale the departure of GDD from GDD crit with GDD length , and to have a phenology factor, φ GDD , that ranges from 0 to 1: When φ GDD = 1, then the Phenostatus cd switches to 3, peak foliage.Full or peak foliage may also occur when carbon allocation to foliage reaches the maximum supported by the available sapwood.Fall senescence (Phenostatus cd = 4) can commence in response to shortening day length ("fall") and decreased air temperature, in a modification of White et al. (1997) and Jolly et al. (2005).Leaves start dropping once air temperature or day length decreases down to threshold values (i.e., T max and ld max ) proportionally to the temperature decreases in the fall as in Eq. ( 5).Full senescence finally occurs when air temperature or day length decrease further down to the minimum thresholds (i.e., T min and ld min ).The phenological factor with respect to air temperature, φ T , is T max , T min , ld max and ld min are constants, with values provided with references in Table 2.

Cold deciduous herbaceous plants
The phenological status of cold deciduous (annual or perennial) herbaceous plants is well captured with functions based on soil temperature (TS), while that of cold deciduous woody plants with air temperature (White et al., 1997).Similarly to Eqs. ( 1) and ( 4) for cold deciduous trees, the soil growing degree days (SGDD) of soil temperature (TS 10 ) are calculated with the base temperature constant (TS base ) of 0 • C. Grasses generate leaves once SGDD exceeds its PFT-dependent critical number (SGDD crit ) and the phenology factor for SGDD, φ SGDD , becomes 1 or greater, as follows: While White et al. (1997) derived SGDD crit as a logistic function of mean annual soil temperature, here we simplify it with three different numbers for different grass types as provided in Table 2.The parameters for φ SGDD were fit to observations at Barrow, Alaska, for arctic C3 grass; the values for C3 and C4 grasses are drawn from White et al. (1997).
Grasses begin fall senescence in response to decreased soil temperature.Leaves start dropping once soil temperature decreases down to a given threshold, TS max , and grasses complete senescence when soil temperature decreases further down to the critical threshold, TS min : β = 1 when the plant is unstressed, and β = 0 at the wilting point.For six soil layers in the LSM, β is calculated for the soil moisture in each layer and averaged weighted by layer thickness and relative root mass fraction, giving the overall β experienced by the plant.
The phenological factor for water stress, ϕ β , is determined by a linear response to the 10-day running average (Foley et al., 1996) of water stress, β 10 , to β max and β min , which represent similarly 10-day running averages of water stress experienced before the onset of drought-induced senescence and at full senescence: When β 10 goes below a minimum (β min ), plants completely senesce in response to drought (ϕ = 0); when β 10 is above a maximum (β max ), plants do not experience drought (ϕ = 1); when β 10 is between β min and β max , the sensitivity of plants to water availability is controlled by the resistance factor (β resis ).The values of s * , s wilt , β min , β max and β resis are provided in Table 2.The leaf-on cue for drought deciduous trees is the same as that for cold deciduous trees, while for grasses the cue is sufficient soil moisture.

Frost hardening in evergreen cold-climate plants
Boreal plants undergo winter frost hardening, which involves physiological changes to protect the plant from chilling injury and freezing injury.Following Repo et al. (1990), the state of frost hardiness S h ( • C) is modeled as follows: where τ is a PFT-specific time constant and the term a•T 10 +b is the stationary frost hardiness, where a and b are PFTspecific parameters for the linear relationship between stationary frost hardiness and air temperature (Hanninen and Kramer, 2007).S h can be thought of as an aggregated measure of the state of the physiological leaf processes that determine the photosynthetic capacity (Mäkelä et al., 2004).
The state of frost hardiness is then used to adjust the maximum photosynthetic capacity V c,max , which is an approach similar to the work of Mäkelä et al. (2006).However, we convert from S h to a dimensionless factor that can take values from 0 to 1.This frost hardiness factor β frost is expressed as where T 0 is a threshold value of cumulative mean temperature at which photosynthesis starts and S h,max is the maximum value of S h (see Table 2 for constants).We implement the first-order Euler scheme to solve Eq. ( 10) and the resulting β frost is used to adjust V c,max .

Fluxnet sites
The Ent TBM was evaluated at five Fluxnet sites, including Morgan-Monroe State Forest (MMSF), Harvard Forest, the Vaira Ranch, the Tonzi Ranch and Hyytiäla, as briefly mentioned above (Table 3).From all sites, data from the flux tower systems were available.Meteorological driver data include radiation, precipitation, air temperature, air pressure, humidity and wind used to drive the model.Soil moisture and temperature measurements were also used to drive the Entstandalone simulations.Flux data includes net ecosystem exchange (NEE) and evapotranspiration (ET), which were used to evaluate the simulation results.Among sites, data availability, such as LAI, varied and suited different types of model simulations as described in detail in the next section.
The MMSF, located in Indiana, USA (Schmid et al., 2000) (latitude: 39.32315 • , longitude: −86.413139 • ) is an extensive, managed temperate broadleaf deciduous forest with a total area of 95.3 km 2 .The area is covered primarily by a secondary successional broadleaf forest within the maple-beech to oak-hickory transition zone of the eastern deciduous forest, dominated by sugar maple and tulip poplar.LAI measurements at 5-14 day intervals during the growing season were available for 1998-2001.
Harvard Forest (latitude: 42.5313 • , longitude: −72.1898 • ) is an eastern temperate mixed forest dominated by deciduous trees.The area surrounding the flux tower is dominated by red oak and red maple, with scattered stands of eastern hemlock, white pine and red pine.About one-third of the existing red oaks were established prior to 1895, another onethird prior to 1930, and the rest before 1940, so the stand is 75-110 years old (Urbanski et al., 2007).O' Keefe (2000) provides the leaf phenology of Harvard Forest.The timing of spring's leaf development and fall's leaf fall have been recorded for permanently tagged individuals in the field since 1991.The leaf development and senescence data in percent of final leaf size have been used to obtain "observational" LAI based on the maximum LAI in the model, The sites are less than 3 km apart.The grasses of both sites are C3 annual species whose growing season is during the winter to spring wet periods.Deciduous blue oaks dominate the savanna overstory of the Tonzi, with a growing season overlapping the grasses during the spring and continuing through the summer drought.In these sites, LAI measurements were available along the tower footprint for 2001 (Kiang, 2002).Hyytiäla (latitude: 61.8474150 • , longitude: 24.294770 • ) in Finland is situated in needleleaf evergreen forest dominated by Pinus sylvestris (Scots pine), in which the phenological behavior of interest is frost hardening.The climate is boreal.Flux measurements and soil moisture and temperature are available.For seasonal LAI, we used the site investigator's description of a constant minimum all-sided needleleaf LAI (75 % of maximum) in January-May, linear increase over June to its maximum of 3.9, remaining at the maximum LAI during July-September, linear decline to its minimum in October and a constant minimum LAI in November-December (Pasi Kolari, personal communication, 2007).

Experiment design
We performed a series of numerical experiments with Ent in different model modes in order to evaluate leaf seasonal dynamics, including leaf phenology and related water and carbon fluxes.We performed simulations for each site with observed soil moisture (hereafter denoted "Ent" mode), and LSM-modeled soil moisture ("LSM" mode); and with observed LAI (without allocation of assimilated carbon to growth; "oveg") and dynamically modeled LAI (via carbon allocation; "dveg"), giving four experiments: Ent-oveg, Entdveg, LSM-oveg and LSM-dveg (Table 4).In the biophysicsonly mode, the observed LAI is prescribed and related active carbon allocations are calculated according to that LAI.In the "active biomass" phenology mode, the leaf phenology and active carbon allocation are dynamically simulated.For MMSF and Harvard Forest, the model was forced with 6 and 9 years' worth of drivers, respectively.In these two sites, continuous soil moisture measurements throughout the rooting depth were not available, so only Ent-LSM simulations were performed.For Vaira, Tonzi and Hyytiäla, the model was forced with a year's worth of tower-measured meteorological drivers as well as observed soil temperature and moisture.
For the Ent versus LSM simulations for annual grass phenology, it was necessary to tailor the soil moisture stress parameters to the different metrics of soil moisture.The phenological timings of grasses depend on the soil moisture condition while an LSM-derived soil moisture is a modelspecific index of soil wetness, not a physical quantity that can be directly validated with field measurements (Koster et al., 2009).The thresholds for the root water stress factor (β in Eq. 8) that was used to model drought deciduous behavior of grasses (volumetric soil moisture at onset of stress and at wilting point) were derived from the observed soil moisture and fluxes, such that these parameters were in a sense tuned to the site as well as to the type of soil moisture measurement.In this study, we therefore tuned the parameters for LSM to better capture the phenological behaviors.
For diagnostics for model performance, we examined observed monthly LAI and monthly sums of GPP, ecosystem respiration (R E ), net ecosystem productivity (NEP = GPP − R E ) and total ET.For potentially waterlimited sites, we examined the modeled volumetric soil moisture and Ent's plant water stress factor.For observed R E , the values are inferred from nighttime respiration and its sen-sitivity to soil temperature, while the modeled values result from both autotrophic and soil respiration.Soil carbon as a driver of soil respiration was initialized from site-measured soil carbon, with litterfall from the model as input on a daily basis (soil carbon was not driven to equilibrium).

Phenology
We evaluated the model performance for cold deciduous woody plants at two sites: MMSF in Indiana and Harvard Forest.
Figure 3 and Table 5 show the simulated variations of the phenological factor (ratio of LAI to the maximum LAI of the year) in comparison to observations.First, it is clear that the gradual nature of changes in LAI during spring and fall were not captured in the model.The phenological factor serves as an on/off cue between environmental thresholds, while growth rate with the ED scheme is limited only by carbon availability, for which reserve carbon is generally not limiting in trees (Sala et al., 2012) or in grass seeds (William Parton, personal communication, 2008).At both sites, the inter-  annual variations of leaf-on timings in the spring were better captured than those of the leaf-off timings in the fall.At Harvard Forest, the dates with the elongation factor of 0.5 in spring showed a correlation coefficient (R) of 0.85 and a root mean squared error (RMSE) of 3.00 days, while the dates with an elongation factor of 0.5 in fall showed R of 0.04 and a RMSE of 15.09 days.

Fluxes
In MMSF, the predicted NEP reasonably followed the observed NEP (Schmid et al., 2000;Dragoni et al., 2007) with correlation coefficients ranging from 0.86 to 0.94, while the peak NEP in summer was slightly underestimated compared to the observed (Fig. 4, Table 6).However, both GPP and R E were more extreme in the model compared to the Fluxnet data product.
In Harvard Forest, the default simulations (LSM-dveg and LSM-oveg) showed underestimated NEP compared to the flux tower observations due to simulated water stress (Fig. 5).As it is known that the cold deciduous plants in Harvard Forest do not experience water stress, no root water stress (β = 1 in Eq. 8) is assumed for additional simulations (LSM-dvegNS and LSM-ovegNS).With the prescribed water stress factor of 1, the model captured the observed NEP reasonably and overestimated GPP and RE compared to observations, similar to MMSF simulations.
The ET values in both LSM simulations were overestimated compared to the flux tower observations in MMSF  and Harvard Forest.These discrepancies might be attributed to both model and data errors.In the model, the higher estimated GPP (although we cannot confirm this) may lead to the overestimated ET to some extent, since higher photosynthesis corresponds to higher canopy conductance and hence more transpiration.In addition, it is well known that eddy flux measurements do not close the energy balance (Wilson et al., 2002).The sum of latent, sensible and ground heat is generally smaller than the net shortwave radiation, which is often caused by measurement errors of latent heat (i.e., ET) and sensible heat (Aranibar et al., 2006), leading to imbalance in measured net radiation of as much as 20 %.The LSM-simulated peak ET is within approximately 70 % of measurements.

Phenology
We evaluated the model performance for drought deciduous herbaceous and woody plants at two sites, the Vaira Ranch and Tonzi Ranch in California.As shown in Fig. 6, at both sites, the timings of C3 annual grasses' green-up and senescence are mainly controlled by soil moisture in a Mediterranean climate, in which precipitation and temperature are seasonally out of phase.Grasses are active during the winter rains but slightly cold-limited in activity and then, with spring warming, growth and activity increase, followed by rapid senescence that closely tracks soil moisture dry down in the late spring and full senescence by the beginning of the dry, hot summer.At the Tonzi Ranch, the oaks have the op- posite seasonality to the C3 grasses.The oaks leaf-out at the end of winter rains around March, when grasses have reached their peak, and then the trees start gradually losing their leaves around the beginning of July due to drought stress.Their complete leaf-off appears to be cued by the November cold or fog, but this latter cue would not be considered a stress factor and is not well understood.At both the Vaira and Tonzi ranches, Ent-dveg and LSM-dveg reasonably captured these phenological timings (Fig. 6).The growth rate for herbaceous plants (i.e., increase in LAI during the growing season) reflected the net carbon assimilation for each day and slightly lagged observations at the beginning of the growing season in the model.Simulated soil moisture clearly decreased much more slowly in LSM-dveg during the late spring dry down compared to the observed volumetric soil moisture that was used to drive Entdveg.

Fluxes
For carbon fluxes at the Vaira Ranch, the model simulations generally followed the observed seasonality, although the late leaf-off in LSM-dveg leads to significant overestimation of carbon uptake, and the observed abrupt increase in R E in the beginning of the growing season was not captured in all cases (Fig. 7, Table 6).Xu and Baldocchi (2004) suggest that the large pulse of R E is the consequence of quickly stimulated microbial activity in decomposition after rain events during the dry season.In the Ent TBM, the soil moisture dependency of decomposition is parameterized as a linear function of soil saturation percent (S) with a plateau when S > S opt (70 %).This response is derived from raw data of soil respiration responses to temperature and moisture in grassland and winter wheat soils from Del Grosso et al. (2005).Most likely, the damped response is because the Ent TBM does not model a separate litter layer on top of the soil and because litter quality may not be well parameterized to allow for fast turnover.As this is a soil model issue, further analysis is worthy of a separate study.
At the Tonzi Ranch, the simulated NEP resulted in a RMSE of ∼ 0.4 µmol m −2 s −1 , as compared to the observed flux (in Fig. 7 and Table 6).During the late spring soil moisture dry-down period, the grasses senesced and the oaks retained their leaves.The oaks started reducing their carbon assimilation due to water stress, as the Ball-Berry slope (m; slope for stomatal conductance) is scaled linearly with the water stress in the model.In reality, the oaks at Tonzi adjust their osmotic potential to maintain their water potential, so their leaf water potential is not linear with soil moisture (Kiang, 2002).Therefore, even with the reasonable LAIs in Ent-oveg, Ent-dveg and LSM-dveg, the underestimated NEP and GPP in the summer are to be expected, lacking a nonlinear response function.Meanwhile, the overestimated LAI in LSM-dveg clearly led to overestimated NEP and GPP.Furthermore, we found the soil biogeochemistry model did not capture the soil respiration pulses after the rainfall, as in Vaira.
The model reasonably captured the observed seasonality of ET with an R of ∼ 0.9 in Vaira and ∼ 0.8 in Tonzi, while the R values for carbon fluxes were much lower.The water fluxes were not much different between LSM-dveg vs. LSMoveg, while the carbon fluxes were significantly different due to different LAIs between the two.The differences in transpiration, resulting from different LAIs, were compensated by evaporation, leading to a relatively small discrepancy in ET between the two experiments.Furthermore, the amplitudes (difference between the maximum and the minimum) of ET were clearly damped in the model, with underestimated peak fluxes during the growing season and overestimated fluxes during the non-growing season.In particular, the noticeable amount of ET occurred during the non-growing season in Vaira, suggesting the partitioning of ET into evaporation and transpiration should be further investigated.

Phenology
At Hyytiäla, the phenological behavior of interest is frost hardening, which lowers photosynthetic capacity in the winter.In comparison to observed LAI, assumed according to Pasi Kolari (personal communication, 2007) and explained in Sect.3.2, simulated LAIs (Ent-dveg and LSM-dveg) (Fig. 8) were almost constant at 4 m 2 m −2 throughout the year, without much decrease during the winter.For evergreen plants, LAI variations in the model reflect the change in foliage carbon balance, as the phenological factor for evergreens remains 1 all the time.Thus, the relatively constant LAIs mean no significant carbon losses during the winter in the model.Based on additional Ent-dveg and LSM-dveg without frost hardening (not shown), we found that such discrepancy in LAI between observation and simulation itself did not influence the predicted water and carbon fluxes noticeably.

Fluxes
Modeled frost hardening in the spring improved the predicted seasonality of NEP markedly in both Ent and LSM simulations (Fig. 8, Table 6).Frost hardening suppressed photosynthetic capacity during the winter (particularly in February-April) and therefore GPP and NEP.It also suppressed transpiration and thus ET, but a relatively small difference in ET was detected between the simulations with and without the frost-hardening scheme as the RMSEs with observations were 7.88 and 7.89 mm s −1 , respectively (Table 6).
With regard to the differences between the Ent-standalone and Ent-LSM models (Ent-dveg vs. LSM-dveg), we found the magnitude of NEP was overestimated in Ent-dveg due to high simulated GPP and underestimated in LSM-dveg due to low soil moisture.During the growing season, the observed volumetric soil moisture was above ∼ 0.35 m 3 m −3 , and the resulting root water stress factor was 1 (completely unstressed) most of the time in Ent-dveg (driven with the observed soil moisture and temperature).However, the predicted volumetric soil moisture was below ∼ 0.25 m 3 m −3 during the growing season in the top three soil layers and the plants roots experienced an average water stress factor of 0.68.Such underestimated soil moisture in the Ent-LSM led to low estimates of NEP.

Discussion
Our experiments show that phenological timing of leaf-out and senescence can be fairly well captured within 10 days or better of observations for deciduous or annual vegetation when based on cumulative weather statistics (e.g., air and soil temperature, growing degree days and day length) derived from observations in the literature.However, the response to soil moisture is sensitive to whether deep root water access is represented to offset soil moisture stress in shallower soil.Also, the soil moisture response must be tuned to the given measure or land model, because soil water content as simulated at the spatial resolution of a land surface hydrology model does not correspond well with any field measure of soil moisture (e.g., volumetric water content, matric potential and pre-dawn water potentials).Stomatal conductance and soil respiration are sensitive to soil moisture stress and hence subject to inaccuracy dependent on the soil moisture representation.Meanwhile, we uncovered weaknesses in the rep- resentation of particular vegetation processes -autotrophic respiration and ED-based carbon allocation -that, besides differences in simulated LAI at one site, were the primary causes of differences from observed NEP.

Drought deciduousness
In Vaira grassland and Tonzi savanna, the phenology parameters which are based on the plant water stress factor (a function of soil moisture), were derived from the site observations of volumetric soil water content (Eq.8) and perform well with observed soil moisture in Ent but not with simulated soil moisture in the LSM.The GISS LSM model predicted the same seasonal trends of soil moisture but higher in magnitude and lower in variability than the observations.Koster et al. (2009) point out that simulated soil moisture is a modelspecific quantity and thus can be considered as an "index" of the moisture state.The specific evaporation and runoff forin addition to model-specific soil parameters such as porosity, hydraulic conductivity, wilting point and layer depth define a dynamic range of soil moisture simulated by the specific model.Therefore, the true information content of soil moisture data lies not necessarily in their absolute magnitudes but in their time variability.
Therefore, the current approach using the absolute soil moisture value for water-limited phenology parameterization could be improved by properly mapping the soil moisture values from the field sites into those in the model, or by using the surrogates for the soil moisture, such as VPD, as suggested by Jolly et al. (2005).However, Stöckli et al. (2011) note that VPD may not be a good indicator of deep soil moisture.
For the trees at MMSF and Harvard Forest, LSMsimulated water stress where the plants should be unstressed indicates that calculating the water stress factor by weighting by root depth distributions does not accurately reflect how trees actually access water.Deep roots generally supply water when shallow layers are dry, and many trees perform hydraulic lift.A future revision of the Ent water stress scheme will account for the ability of plants to preferentially access soil moisture at any depth in the root zone, such that soil moisture stress is not a simple weighted average through the root profile.
While the Fluxnet data have recently been widely used to evaluate the DGVMs and LSMs, we still find the need for more comprehensive measurements at the sites.Specifically, it was very difficult to have continuous soil moisture and temperature together along with measurements from eddy covariance towers; also, the detailed tree surveys were not always available.

Cold deciduousness
For cold deciduous trees, we used the growing degree days and chilling requirements in spring phenology (Botta et al., 2000) and temperature and photoperiod in fall phenology (White et al., 1997;Jolly et al., 2005).While we have taken a widely used approach, some recent studies suggest other possible approaches.For spring phenology, the importance of the photoperiod has been pointed out in recent studies (e.g., Körner and Basler, 2010;Migliavacca et al., 2012).Körner and Basler (2010) suggested that when the chilling requirement is fulfilled, plants become receptive to photoperiod signals and such sensitivity to the photoperiod is found in late successional species in mature forests.For fall phenology, Delpierre et al. (2009) used chilling a degree-day photoperiod to model leaf coloring change for deciduous trees in France, and Yang et al. (2012) and Archetti et al. (2013) found the model suitable for New England, US, with different parameter fits.In general, despite agreement about overall climate cues for cold deciduousness, further work is needed to uncover site-independent parameterizations.

Photosynthesis and respiration parameters
In this study, site-specific parameters were used according to the data availability.As in Appendix D, some of parameters are generic for PFTs and some are site-specific.For the model to be utilized at the global scale, further exploration is required to determine geographic variation in parameters and possible climatology-based prediction of parameters.Model parameters for biophysics or ecosystem models have been inferred with various mathematical techniques, such as a Monte Carlo simulation (Kleidon and Mooney, 2000), data assimilation with Kalman filtering (Mo et al., 2008;Stöckli et al., 2008), optimization with the Marquardt-Levenberg method (Wang et al., 2007) and optimization with the simulated annealing method (Medvigy et al., 2009;Kim et al., 2012).
In general, vegetation biophysics models can replicate observed canopy fluxes of CO 2 well when the vegetation structure is well specified, but the same net flux can be predicted from different levels of gross assimilation versus respiration.The main biophysical parameters common to most models are the maximum leaf photosynthetic carboxylation rate, V cmax ; autotrophic respiration as a function of biomass, temperature and activity; and leaf litter quality, such as lignin content, for soil respiration.While V cmax may be precisely measured for a leaf, its value can be highly variable within a plant and seasonally.Currently, in the Ent model, V cmax is only variably with PFT and temperature, and the intrinsic quantum efficiency for J max , the maximum leaf photosynthetic electronic transport, is constant.The seasonal variation of V cmax , J max and SLA (specific leaf area) could be introduced, pending better mechanistic understanding.A simple approach based on photoperiod such as in Bauerle et al. (2012) would be possible.
Autotrophic respiration can range ∼ 30-80 % of annual GPP for different plant types (Falge et al., 2002).These parameters, however, may not extrapolate to the global scale, so future study is necessary to investigate global variation in parameterizations.In general, respiration is poorly understood and cannot be modeled fully mechanistically but must rely on bulk parameterizations that effectively integrate numerous processes.Researchers have attempted various approaches to grouping some respiratory fluxes (Amthor, 2000;Cannell and Thornley, 2000) as responsive to different drivers, though there is as yet no generally accepted scheme.In Ent, the streams are maintenance respiration, which is a function of biomass and responsive to temperature, "light growth respiration" from photosynthetic activity and "biosynthesis respiration" from growth or turnover of plant tissues.
In Ent, using site-specific parameters for leaf photosynthetic capacity, V cmax , constant throughout the canopy, we observed a tendency toward higher GPP and higher ecosystem respiration, R E , compared to that inferred from tower observations when nighttime respiration temperature response is used to estimate R E .These extremes in the two components of the net flux are not necessarily unreasonable, since the Fluxnet respiration product could be underestimated.The R E data products we used were modeled, as usual, with an exponential equation to fit the measured nighttime CO 2 flux as a function of soil temperature (Schmid et al., 2000).Such an estimate excludes daytime root respiration, which increases with photosynthetic activity (Tang and Baldocchi, 2005;Tang et al., 2005).With regard to GPP, recent oxygen isotope work suggests that global gross primary productivity is higher than traditional estimates (Welp et al., 2011).It is a well-known problem in ecosystem science that GPP and respiration cannot be directly partitioned through current measurement methods for net ecosystem exchange, although there are hopes for a solution now possible with measurements of solar-induced fluorescence (van der Tol et al., 2014).

Carbon allocation/growth scheme
We encountered deficiencies in the carbon allocation/growth scheme that we adopted from the ED model.Although the current carbon allocation and growth scheme results in LAI that is reasonable, with some phenological timing issues as noted, the maximum LAI is achieved thanks to a cap on LAI by allometric relations to stem structure and plant density, while the rest of the plant carbon balance is not realistic, particularly with regard to rate of LAI growth, amount of seasonal sapwood growth and conversion to heartwood, accumulation of carbon reserves and allocation to reproduction.The on/off cues of the Ent phenological factor for cold deciduous trees results in an unrealistic fast full leaf-out, which could be rectified by introduction of a physically based cell growth elongation factor (Lockhart, 1965).We also found it would be more realistic to make carbon allocation to each live pool independent.The ED scheme allocation to one live biomass total and then partition among the live pools can lead to unrealistic behaviors for sapwood patterns during spring growth and fall senescence, due to a partitioning scheme for live carbon that does not account for the different seasonal behaviors of each live pool.Finally, reproduction in ED is currently a fixed fraction of assimilated carbon, which is problematic in the plant's overall carbon balance as a large sink.Recent studies show that reproduction relies heavily on stored carbon, which often accumulates over more than a year, such that growth of other plant tissue is never carbon limited while large stores are kept in reserve.The ED scheme relies on the plant using nearly all stored carbon for deciduous plants each year.Introducing reproductive allocation based on thresholds proposed by Sala et al. (2012) would help rectify Ent's simulated plant carbon balances such that trees are not always reaching the limit of carbon starvation.Besides respiration, plant carbon allocation is currently still poorly understood.However, recent studies with carbon tracers (Epron et al., 2012a, b) are yielding new insights that could be used to improve growth schemes that continue to be a weakness in dynamic global vegetation models.

Conclusions
In this study, we evaluated the Ent TBM focusing on the seasonal dynamics of vegetation leaf as well as carbon and water fluxes.In particular, we took a process-based approach, evaluating the Ent-standalone model with observed LAI and Ent's prognostic active growth submodel with observed soil moisture as well as coupled to the LSM model for prognostic soil moisture, allowing us to identify parameterizations that need to be improved.For herbaceous PFTs whose phenological timings depend on soil water availability, it is inevitable to find errors in phenological timing in Ent-LSM simulations due to the discrepancy in simulated soil moisture in the LSM.Also, the predicted LAI of herbaceous PFTs in Ent directly reflects the amount of assimilated carbon on the day and vice versa as herbaceous PFTs allocate assimilated carbon only to active compartments (as they have no structural tissue) and thus any errors in phenological timings propagate into errors in biophysical processes.For tree PFTs, the Ent soil moisture stress scheme should be improved to allow for deep soil moisture access to override stress that might result from weighting shallower dry soil layers too strongly.
This study evaluated the phenology and resulting seasonality of fluxes in the limited number of sites, including four different PFTs.The Ent PFTs not tested in this study include deciduous needleleaf plants, evergreen broadleaf plants, shrubs, arctic grasses and crops.Future work will involve determining the efficacy of these PFT parameterizations at the global scale and the possibility of developing some of these parameters as functions of local climate as obtained from either reanalysis data or from GCM climatology.In addition, we have identified deficiencies in the carbon allocation scheme from the ED model that can be rectified in future revision of Ent's growth submodel.Future work will include development of phenology and allometry parameter sets that are robust at the global scale and soil moisture stress accounting for deeper soil access.In addition, due to how ED allocates biomass to all live pools (e.g., foliage, sapwood and fine roots) combined, rather than allowing for separate dynamics, alternative carbon allocation schemes that partition the dynamics of the live tissues must be developed for realistic plant carbon balances.
This work sets the foundations for coupled land carbon-GCM simulations that can utilize height-structured canopy data from remotely sensed lidar to reduce uncertainty in pre-dictions of the land carbon balance through tighter links between seasonal growth dynamics geometrical and biomass allometry of vegetation canopies.Because the model at the global scale will involve a community of users that will continue to identify parameter sets applicable for more climatically diverse distributions of the Ent TBM's PFTs, this paper is also written to serve as a detailed reference for these users, to achieve appropriate interpretation of model results and parameter adjustment.1.1 1.1 0.25 0.25 0.25 0.25 0.25 0.25 0.8 0.8 0.9 0.9 0.9 0.9 0.9 0.25 b 0.4 0.4 2.0 2.0 2.0 2.0 2.0 2.0 0.4 0.4 0.9 0.9 0.9 0.9 0.9 2.0 The soil biogeochemistry submodel of Ent utilizes a slightly modified version of the CASA' biosphere submodel originally implemented in the NCAR (National Center for Atmospheric Research) LSM and CSM 1.4 (Bonan, 1996;Randerson et al., 1997;Fung et al., 2005;Doney et al., 2006), which itself is a modified version of the original NASA CASA biosphere model (Potter et al., 1993).The soil model determines terrestrial soil carbon pools and CO 2 fluxes from microbial respiration.

B1 CASA structure
The soil biogeochemistry model consists of three litter C and N pools and nine soil C and N pools, as in CASA'.The pools are currently only simulated for the top 30 cm soil depth.This layer accounts for nearly all observable soil respiration fluxes to the atmosphere but not for full long-term carbon stocks in deeper soil.Simulating soil carbon down to 100 cm and deeper would allow comparison to existing global data sets of soil carbon and root depths (Batjes, 1996a, b;Jackson et al., 1996).Figure B1 shows these 12 pools.Ent has an optional 30-100 cm deep soil layer that is not run in the current paper.
The various pools currently have fixed C : N ratios and turnover times, listed in Table B1.The pools gain carbon and nitrogen from transfers from other pools and losses to respiration and transfers to other pools.These transfer and respiration fractions are listed in Table B2.
Soil micrometeorological conditions for the soil layers must be extrapolated from the soil layering scheme of the land surface model.For example the GISS land surface hydrology has a six-layer soil scheme with geometrically increasing layer thicknesses with depth (Rosenzweig and Abramopoulos, 1997), so soil temperature and moisture for the soil biogeochemistry layers are calculated through a weighted sum for the upper 30 cm.
In addition to the transfer coefficients in Table B2, three other rate coefficients are used (following Randerson et al., 1997): These are simply decomposition rate adjustment factors for soil microbial, slow and passive pools (respectively) for crops only; their value for all other PFTs is 1.

B2 Soil module interface with vegetation
Physical inputs to the soil module from the land surface hydrology are volumetric soil moisture, soil temperature and soil texture (percentage of clay, sand and silt).Biological inputs consist of leaf, root and wood litter (Fig. B1).Model  a From the original CASA code (Potter et al., 1993).b From the CASA' code (Doney et al., 2006).
outputs are soil C (and N, not used) pools and soil CO 2 flux.Ent calculates litterfall carbon from the leaf area times the specific leaf area.The relevant PFT-dependent litter parameters (leaf, fine root and wood turnover times, litter C : N ratios, specific leaf area and lignin contents) from Ent are listed in Table B2.In addition to these parameters, a parameter representing the inverse of the residence times of the litter pools, denoted annk lit (in units of year −1 ), was calculated as the inverse of lrage (leaf and root litter age) for leaf and root litter or of woodage (stem litter age) for wood litter (Potter et al., 1993).

B3 Temperature and moisture responses of soil respiration
We replaced the CASA' temperature and soil moisture responses of soil respiration with new functions derived from new fits to field data collected by Del Grosso et al. (2005).The Ent TBM temperature response of soil respiration is a simple piecewise linear model that increases up to 30 • C and then flattens.In reality, the response to temperature is exponential up to a certain optimum then declines, but a linear representation was chosen because it reduces the computational time compared to that required for calling an exponential function, and tests on field data show adequate performance for the purpose of predicting respiration fluxes and soil carbon pools (unpublished).At high soil temperatures, soil moisture stress usually also occurs but, because no measurement data were available for respiration at temperatures above 30 • C, the Ent model response does not represent a decline in soil respiration at high temperature.The linear tem- perature response of soil respiration is More realistically, the temperature response is in nature an exponential response, so if there are no computational constraints the following Q 10 function as formulated in the original CASA' should be used: where Q 10 has a typical value of 2.0.The Ent TBM moisture response of soil respiration is similarly a piecewise linear model that rises from 0 at zero soil moisture to 1.0 at a relative extractable water content (REW) of 0.7, where REW is the fraction of saturation above the hygroscopic point.Because there are no good functions for calculating the hygroscopic point based on soil texture, we estimate the hygroscopic point as half of the wilting point.We note that it would be more precise to model the soil moisture response as an optimality curve that rises from the soil hygroscopic point (minimum for microbes) rather than wilting point (for plants) to some optimum and then declines as pore space becomes saturated and obstructs the flux of gases.However, because of a lack of good algorithms to calculate the soil hygroscopic point for different soil textures, we use this version of Ent relying on the wilting point as the point of minimum available soil moisture.We may later introduce a simple linear decline of the soil moisture response with saturation; however, at present we have no data on the response The linear soil respiration temperature and moisture response functions are plotted in Appendix Figs.B1 and B2, along with the original CASA' responses and those of Del Grosso et al. (2005), whose data were re-analyzed to generate the Ent response functions.

Appendix C: Allocation
The labile carbon reserves in Ent are allocated into different plant biomass pools, including foliage, sapwood, heartwood, fine root and coarse root.In addition, turnover of tissue due to background litter fluxes is replenished from the carbon reserve pool.In nature, plants may allocate biomass to different compartments in response to many different controlling factors, such as light availability and water availability, which alter, for example, root : shoot ratios.Among various carbon allocation modeling approaches with different complexities, many DGVMs take a simple approach to model carbon allocation via empirical and allometric relationships, a traditional "static allocation" approach (Foley et al., 1996;Sitch et al., 2003) while some models parameterize the dependency of carbon allocation on resource availability, "dynamic allocation" approach (Friedlingstein et al., 1999;Arora and Boer, 2005).Although carbon allocation varies with plant characteristics such as size and age and environmental conditions, the static allocation approach may be justified for models operating at large scale.If plant productivity is assumed in a steady state, carbon allocation is likely to be in a steady state.Also, spatial variability in environmental factors and their effects on allocation can be averaged.However, the fixed allocation approach is limited in long-term simulations as it lacks response to environmental changes such as climate change and elevated atmospheric CO 2 (Franklin et al., 2012).However, recent models of dynamic allocation have been difficult to constrain due to a dearth of observations.In the Ent TBM, the allocation submodel takes a traditional approach of static allocation, based on allometric relationships between different pools.Modified from approaches of the ED models (Moorcroft et al., 2001;Medvigy et al., 2009), the scheme allocates the labile carbon to different biomass pools according to empirical allometric relationships and leaf phenological status on a daily basis.

C1 Active biomass
The biomass within each plant is partitioned between an active carbon pool and a structural carbon pool.The active biomass pool (B active ) (kg biomass /cohort) is subdivided into foliage (B fol ), sapwood (B sw ) and fine roots (B froot ) which turn over at different rates, while the structural pool (B structural ) consists of heartwood (B hw ) and coarse roots (B croot ).Grasses do not have the structural pool.The labile biomass (B lab ) assimilated on the same day is allocated to the active carbon pool to maintain the size of foliage, sapwood and fine root tissues given their turnover rates and to accumulate the active carbon up to its maximum.
Thus, the time change of the active pool can be written as where B max active is the maximum active carbon of each plant, which is determined according to the maximum foliage carbon according to the size of the plant, CB d is the daily plant carbon balance (i.e., sum of NPP on one day).Then, the allometric relationships are used to subdivide the active biomass into its components.The foliage biomass is determined according to its phenological status (ϕ), ranging from 0 (for full senescence) to 1 (full leaf-out) as a proportion of full-leaved foliage biomass, B * fol , so that B fol = ϕB * fol .Both the fine root and sapwood biomass are also determined according to their proportional relationships to B * fol .A constant empirical proportionality for fine root (q fr ), assumed to be 1, is as follows: The sapwood biomass is determined according to the pipemodel theory (Shinozaki et al., 1964), which suggests that the total foliage area is proportional to the sapwood crosssectional area.The ratio between full-leaved foliage area and sapwood area is assumed to be 3900 (m 2 foliage m −2 sapwood ).This value is adopted from the value used in ED1 (Moorcroft et al., 2001), which follows Rending and Taylor (1989), giving the ratios of foliage area to sapwood area a range from 3900 to 14 000.These assumptions result in the following relationship: where ρ sw is the sapwood density (kgC m −3 sapwood ) and SLA is the specific leaf area (m 2 foliage kgC −1 ) for each PFT, provided in Table 1.ρ sw is taken to be 500 (kgC m −3 sapwood ) (i.e., 0.5 kgC kg −1 biomass × 1000 kg biomass m −3 sapwood for very hard wood).However, we note that there are departures from these constant values.The fraction of dry biomass that is carbon in spruce wood is typically 0.48 (Payne, 2002).Also, Schneider et al. (2011) find the foliage to sapwood area ratio to be closer to 500-600 for Jack pine, with higher values toward the interior of the sapwood that serves older foliage.Calvo-Alvarado et al. (2008) find an increasingly linear relationship between height and foliage area/sapwood area for Costa Rican rainforest trees, ranging from 500 to 1500.A consistent rule for this variation has yet to be identified, but it may vary with wood density and anatomy.
Finally, B fol is related to LAI (m 2 foliage m −2 ground ), measuring the total leaf (i.e., foliage) area per the projected ground area as where nplant is the population density of cohorts (# plants m −2 ground ) and 0.5 (kgC kg −1 biomass ) is to convert SLA (from m 2 foliage kgC −1 to m 2 foliage kg −1 biomass ).

C2 Structural and reproductive biomass
Growth of structural tissue is handled as follows.If the stored labile biomass is non-zero, the size of the structural pool of woody plants increases according to the empirical allometric relationships and consequently the size of the active pool increases.Here, the partitioning between B active and B structural is written as where DBH is the diameter at breast height and q structural is the ratio of structural growth to active growth.The derivatives are derived from allometric relationships according to plant size (i.e., DBH and height) for woody plants.Note herbaceous plants do not have the structural pool, meaning that DBH = 0, q structural = 0, B structural = 0 and q sw = 0. Also, the plant devotes a fixed fraction (q repro ) of daily carbon to the reproductive pool and the rest to growth of the active and structural pools.q repro is assumed to be 0.3 for woody plants and 1.0 for herbaceous plants, following the assumptions of ED1 (Moorcroft et al., 2001).

Appendix D: Biophysics, allocation, and phenology parameters
See Tables D1 and D2.

Figure 1 .
Figure 1.Schematic diagram of the Ent model.

Figure 2 .
Figure 2. Ent's individual plant biomass pools and geometry.Herbaceous plants exclude woody tissue.

Figure 3 .
Figure 3. Daily simulated (S) and observed (O) phenology: (top) LAI/LAImax, (middle) phenological dates (day of year) for spring leaf-out at percentage of maximum, and (bottom) phenological dates (day of year) for fall senescence in MMS and Ha1.These results show good simulated timing of initial leaf-out and final senescence but lack of the gradual rate of these, such that the maximum leaf-out occurs too soon and the period of peak growth is too long.The gradual behavior could be simulated through a rate constraint.

Figure 6 .
Figure 6.(a) Daily root water stress and (b-c) daily LAI in the Vaira and Tonzi ranches for 2002.

Figure 8 .
Figure 8. Monthly fluxes and daily states in Hyytiäla for 1998: (a) NEP, (b) GPP, (c) RE, (d) ET, (e) LAI, (f) soil temperature and (g) root water stress.Here the observed LAI is assumed based on personal communication with the site investigator, Pasi Kolari (2007).
Figure A1.(a) Cumulative root density profile distributions and (b) probability density distributions in the Ent TBM (modified from Rosenzweig and Abramopoulos, 1997) vs. soil depth increments of the NASA GISS GCM land surface model.

Figure B1 .
Figure B1.Schematic diagram of the soil biogeochemistry submodel of Ent (showing nine soil C pools only; modified from Potteret al., 1993).Surfstr -surface structural pool; Surfmet -surface metabolic pool; Soilmet − soil metabolic pool (fastest to decompose; 20-day turnover time); Soilstr -soil structural pool; Surfmicsurface microbial pool; Soilmic -soil microbial pool; Slow -slowly decomposing pool; Passive -very slowly decomposing pool (500year turnover time).All pools except for the three surf * * * pools are assumed to be present in the two lower soil layers in addition to the top layer.

Figure
Figure B2.(a) Temperature responses of soil respiration in Del Grosso et al. (2005), CASA' and Ent's piecewise linear response; and (b) moisture response of soil respiration in Del Grosso et al. (2005), CASA' and Ent for grassland (Vaira Ranch) soil texture.

Y.
Kim et al.: Variability of phenology and fluxes of water and carbon in Ent TBM to saturated conditions.

a
For all these plant functional types there is a large range of values as well as large variation within a single site and single plant.We therefore have chosen literature values for the Fluxnet sites where available and tuned the value within the literature range for the site.bOleson et al. (2004).c Wilson et al. (2002).dWang et al. (2007).

Table 1 .
Plant functional types in Ent.

3865, 2015 3844 Y. Kim et al.: Variability of phenology and fluxes of water and carbon in Ent TBM 2.4.4 Drought deciduous woody and herbaceous plants
(Rodriguez-Iturbe et al., 2001): available soil water for the plant.In the model, it is determined based on a 10-day running average of the physical time step (∼ half-hourly) plant water stress factor β.This factor is the same used to scale stomatal conductance for water stress and is determined by a linear response between PFT-dependent critical relative soil moisture (volumetric soil moisture/saturated volume) points for the plant at which water stress begins, s * , and at which wilting occurs, s wilt ,(Rodriguez-Iturbe et al., 2001): SeeTable 2 for constant values of TS max and TS min .www.geosci-model-dev.net/8/3837/2015/Geosci.Model Dev., 8, 3837-
i.e., observed LAI = observed percent of leaf development or fall × modeled maximum LAI.The Vaira Ranch (latitude: 38.4066667 • , longitude: −120.950733• ) and Tonzi Ranch (latitude: 38.4316 • , longitude: −120.9660• ) in Ione, California, are located in an open grassland ecosystem and an oak/grass savanna ecosystem, respectively, in a Mediterranean climate of cool wet winters and dry hot summers.

Table 4 .
Types of experiments.

Table 5 .
Correlation coefficients and RMSEs of LAI-based phenological dates between simulations and observations.

Table 6 .
Correlation coefficients and RMSEs of hourly and daily fluxes between simulations and observations.

Table A1 .
Plant functional type parameters for root density distributions.

Table B1 .
Values of C pool parameters: C : N ratio of all 12 C pools (used only to calculate N pools); annk soil -inverse of turnover times of all 9 soil C pools (year −1 ).
(Doney et al., 2006)sfer efficiencies for the physical time step are from the CASA' code as implemented in the National Center for Atmospheric Research (NCAR) climate model (CSM) 1.4(Doney et al., 2006).

Table D1 .
Biophysics parameters for Fluxnet sites in this study.

Table D2 .
Biogeochemical and phenological parameters for Fluxnet sites in this study.