An improved parameterization of leaf area index (LAI) seasonality in the Canadian Land Surface Scheme (CLASS) and Canadian Terrestrial Ecosystem Model (CTEM) modelling framework

Leaf area index (LAI) and its seasonal dynamics are key determinants of vegetation productivity in nature and as represented in terrestrial biosphere models seeking to understand land surface atmosphere flux dynamics and its response to climate change. Non-structural carbohydrates (NSCs) and their seasonal variability are known to play a crucial role in seasonal variation in leaf phenology and growth and functioning of plants. The carbon stored in NSC pools provides a buffer during times when supply and demand of carbon are asynchronous. An example of this role is illustrated when NSCs from previous years are used to initiate leaf onset at the arrival of favourable weather conditions. In this study, we incorporate NSC pools and associated parameterizations of new processes in the modelling framework of the Canadian Land Surface Scheme-Canadian Terrestrial Ecosystem Model (CLASS–CTEM) with an aim to improve the seasonality of simulated LAI. The performance of these new parameterizations is evaluated by comparing simulated LAI and atmosphere–land CO₂ fluxes to their observation-based estimates, at three sites characterized by broadleaf cold deciduous trees selected from the FLUXNET database. Results show an improvement in leaf onset and offset times with about 2 weeks shift towards earlier times during the year in better agreement with observations. These improvements in simulated LAI help to improve the simulated seasonal cycle of gross primary productivity (GPP) and as a result simulated net ecosystem productivity (NEP) as well.

1 Introduction

Biosphere–atmosphere interactions constitute a complex system which plays an important role in the regulation of the climate. These interactions are important determinants governing the physical and chemical properties of the atmosphere as well as the growth of plants and result in the biosphere and atmosphere behaving as a coupled system (Pilegaard et al., 2003). Understanding this coupled behaviour is a key research priority due not only to the important role that terrestrial ecosystems play in modulating the global carbon cycle but also to the significance of land surface characteristics for local and regional climate through biogeophysical effects (Cox et al., 2000; Prentice et al., 2001; Bonan, 2008; Franklin et al., 2016). This growing recognition of the role of land surface vegetation, and its bidirectional interactions with the climate system, has led to ever increasing complexity of the physical and biogeochemical processes that are incorporated in the land surface components of regional and global climate models (Foley et al., 1996; Sitch et al., 2008; Flato et al., 2013). Process-based land surface schemes and ecophysiological models (e.g., Running et al., 1999; Mäkelä et al., 2000; Friend et al., 2007; IPCC, 2013; Sato et al., 2015) simulate atmosphere–land exchanges of carbon, water, and energy and offer tools for understanding vegetation behaviour for the present climate and for projecting vegetation behaviour for future climate scenarios.

The plant canopy is a locus of physical and biogeochemical processes in an ecosystem. The functional and structural attributes of plant canopies are affected by microclimatic conditions, nutrient dynamics, herbivore activities, and many other factors (Asner et al., 2003). Leaves are the point of contact between plants and atmospheric CO₂; an increase in leaf area potentially enhances the opportunity for carbon uptake, albeit at the cost of a greater demand for water (Norby et al., 2003). The amount of foliage contained in plant canopies is one of the most basic ecological characteristics that integrates the effects of overall environmental conditions. Canopy leaf area serves as the dominant physical control over primary production (photosynthesis), transpiration, energy exchange, and other physiological attributes pertinent to a range of ecosystem processes and is therefore a core element of ecological field and modelling studies (e.g., Knyazikhin et al., 1998; Xavier and Vettorazzi, 2004; Aboelghar et al., 2010; Gonsamo and Chen, 2014; Bao et al., 2014; Savoy and Mackay, 2015).

Leaf area index (LAI, defined as the amount of leaf area (m²) in the canopy per unit ground area (m²)) is a dimensionless quantity and therefore can be assessed across a range of spatial scales, from individual plant, a forest stand, or grassland to large regions and continents. Leaf phenology describes the response of leaves to seasonal and climatic changes including the timing of bud burst, senescence (leaf maturity or browning), and leaf abscission (leaf fall) and has been documented in a wide range of literature (e.g., Kikuzawa, 1995; Myneni et al., 1997; Arora and Boer, 2005; Menzel et al., 2006; Parmesan, 2006; Richardson et al., 2010; Dragoni et al., 2011; Smith and Hall, 2016). Leaf phenology is a function of environmental conditions (in particular, temperature, soil moisture, and day length). The structural and adaptive qualities specific to vegetation type also determine the timing of leaf phenological events. Accurate prediction of recurring vegetation cycles as a function of climate is an important feature that vegetation models are expected to reproduce. The timing of bud burst and leaf senescence determines the length of the growing season, and this affects gross and net primary productivities (GPP and NPP), the annual cycle of LAI, and consequently the energy, water, and carbon fluxes. The seasonal progression of LAI also influences canopy conductance (Blanken and Black, 2004), albedo (Sakai et al., 1997), and through its modulation of sensible heat (H) and latent heat (LE) fluxes (Moore et al., 1996) also surface air temperatures (Levis and Bonan, 2004).

Despite its importance, the representation of LAI in terrestrial biosphere models is considered poor (Richardson et al., 2012). Lack of high-quality long-term observations, the use of prescribed LAI, simplified formulations of underlying biogeochemical processes, and coarse spatial resolution have been mentioned as some of the limitations to accurate representation of LAI (Kucharik et al., 2006). Since canopy seasonality is an important determinant of carbon (C) fluxes, poor representation of the seasonal dynamics of LAI can lead to inaccurate estimation of vegetation productivity and consequently the net atmosphere–land CO₂ flux (Ryu et al., 2008).

Non-structural carbohydrates (NSCs) are the primary products of photosynthesis and a key energy source for plant growth and metabolism. NSCs play a central role in a plant's life processes and its response to the environmental conditions (Kozlowski, 1992; Ögren, 2000; Chatterton et al., 2006; O`Brien et al., 2014; Hartmann and Trumbore, 2016; Sperling et al., 2017). Previous studies have suggested that NSCs are stored in all plant organs (i.e., leaf, branch, root, and stem) at different concentrations that vary seasonally and also inter-annually in response to changes in environmental conditions (e.g., Oberhuber et al., 2011; Bazot et al., 2013; Mei et al., 2015). The amount of NSC and its particular allocation to leaves, stems, and roots are considered ecophysiological traits and are among the range of adaptive strategies that plants use (Li et al., 2001; Poorter and Kitajima, 2007; Wyka et al., 2016). Many factors influence leaf NSC content, including nutrient elements (Zotz and Richter, 2006), temperature (Gough et al., 2010), precipitation (Würth et al., 2005), drought (Rosas et al., 2013), and phenology (Chen et al., 2017). NSC concentrations also vary seasonally and this seasonality has been observed in various plant species (e.g., Zhu et al., 2012; Richardson et al., 2013; Saffell et al., 2014). Despite this extensive research, the size and relative contributions of NSC pools across different tree organs and the movement of NSC amongst the tree organs are not well understood (Friedlingstein et al., 1999; Allen et al., 2005; Sala et al., 2012; Mei et al., 2015).

Plant NSC stores can compensate for a carbon or nitrogen shortage when current demand surpasses supply due to the seasonality of plant growth, stresses, or disturbances. The seasonal dynamics of NSC concentrations have been studied in various plant species (e.g., Zhu et al., 2012; Richardson et al., 2013; Saffell et al., 2014). In deciduous plants, when photosynthesis is constrained by limited leaf area and low temperature in early spring, NSC is mobilized from stem and roots to support respiration and tissue growth, resulting in decreased concentrations of NSC in these storage organs (Hoch et al., 2003; Palacio et al., 2007). During the growing season, storage pools are replenished and NSC concentration increases (Teixeira et al., 2007; Klein et al., 2016). Typically, NSC concentrations in storage organs of the short-lived fast-growing species decrease in springtime after bud flush and then increase during the remainder of the growing season. Correspondingly, the storage organs shift from being a NSC source in the early growing season to becoming a sink in the late growing season, maintaining tree survival after the termination of photosynthate flow from aboveground sources to supply energy for stem and root tissues through the winter (Würth et al., 2005; Gough et al., 2010). During periods of limited photosynthesis, such as winter dormancy or drought stress, trees depend solely on stored NSCs to maintain basic metabolic functions, produce defensive compounds, and retain cell turgor (Sperling et al., 2015). For deciduous species, the NSC storage provides the means to jump-start leaf onset by using a part of NSC stores to push leaves out at the onset of favourable weather conditions (e.g., in spring in the Northern Hemisphere). Representation of NSC pools is therefore an essential step for terrestrial biosphere models to better simulate leaf phenology and seasonal variability in LAI.

Models typically tend not to explicitly represent NSC pools (Le Roux et al., 2001) because of the additional complexity that has to be dealt with unless it is necessary to do so for improved physiological realism of the overall model or representation of a particular physiological process. Those with explicit NSC pools (Dick and Dewar, 1992; Sampson et al., 2001; Ogle and Pacala, 2009; Fisher et al., 2010) may choose to assign no metabolic cost associated with NSC pools and treat it simply as a storage pool from which carbon can be used at a future time (Sala et al., 2012), or they choose to explicitly account for the metabolic carbon costs associated with NSC pools. Existing representations of NSC pools in models have aimed to maintain metabolic function through carbon deficit periods (Sampson et al., 2001), achieve trade-offs between seedling growth and survivorship (Kobe, 1997; Fisher et al., 2010), simulate stress tolerance to ozone (Retzlaff et al., 1996), and represent the response of the respiration ∕ photosynthesis ratio to increase in temperature (Dewar et al., 1999). Models have also used mass balance of NSC pools as a trigger for plant mortality (Dick and Dewar, 1992; Cropper and Gholz, 1993; Sampson et al., 2001; Ogle and Pacala, 2009; Fisher et al., 2010). In terms of model structure, models may have only a single NSC storage pool (Cropper and Gholz, 1993; Kobe, 1997; Schaefer et al., 2008; Ogee et al., 2009) or multiple above- and below-ground pools (Le Dizeìs et al., 1997; Berninger et al., 2000; Allen et al., 2005). Biochemically, models vary in what the NSC pools correspond to. Most models do not distinguish between soluble (e.g., sugars) and insoluble (e.g., starch) NSCs, whereas others may make this distinction (Dick and Dewar, 1992; Le Dizeìs et al., 1997; Dewar et al., 1999; Levy et al., 2000; Geìnard et al., 2008; Schaefer et al., 2008). Dietze et al. (2014) provide a review of existing models which incorporate NSC pools within their framework.

Here, we include a representation of NSC pools and the associated parameterizations in the framework of the Canadian Land Surface Scheme–Canadian Terrestrial Ecosystem Model (CLASS–CTEM). CLASS–CTEM exhibits delayed leaf phenology and we attempt to address this issue. In the original model, the simulated global LAI reaches its maximum in August whereas the observed LAI peaks in July (e.g., see Fig. 11 of Anav et al., 2013). The objective of this study is to improve and assess the performance of CLASS–CTEM simulated leaf phenology for broadleaf cold deciduous trees. Model performance is evaluated against in situ measurements from three sites in the FLUXNET data network (https://fluxnet.ornl.gov/obtain-data, last access: 14 March 2018) which provides tower-based meteorological variables used to drive the model as well as observations of LAI, carbon, and energy fluxes.

Figure 1 Schematic representation of the CTEM model after addition of non-structural carbohydrate pools. The arrows in blue show the new non-structural carbohydrate fluxes as shown in Eqs. (5) and (6).

2 Model, data, and methods

2.1 CLASS–CTEM model

A coupled version of the Canadian Land Surface Scheme (v. 3.6; Verseghy, 2012) and Canadian Terrestrial Ecosystem Model (v. 2.1.1, Melton and Arora, 2014) (CLASS–CTEM) is used here. Slightly older versions of both models are currently implemented in the second-generation Canadian Earth System Model (CanESM2; Arora et al., 2011). While CLASS simulates fluxes of energy, water, and momentum at the land–atmosphere boundary, atmosphere–land fluxes of CO₂ and CH₄ are simulated by CTEM. CLASS operates at a typical time step of 30 min and prognostically simulates the liquid and frozen soil moisture and soil temperature for its multiple soil layers (three layers are employed here with maximum thicknesses of 0.1, 0.25, and 3.75 m); the temperature, thickness, and fractional cover of snow; and the temperature and amount of snow and rain on the vegetation canopy. The permeable depth of the soil column may be smaller than the total depth of the soil layers employed; if a layer spans the permeable depth boundary it is subdivided for hydrological calculations. CLASS distinguishes four plant functional types (PFTs) (needleleaf trees, broadleaf trees, crops, and grasses) as shown in Table 1. CLASS calculates net radiation (R_n) based on prognostically calculated land surface albedo and the skin temperature of the land surface (T_s) as

$$\begin{array}{}\text{(1)}& R_{\mathrm{n}}=\mathrm{SW}\left(\mathrm{1}-\mathit{\alpha }\right)+\mathrm{LW}-\mathit{\sigma }{T}_{\mathrm{s}}^{\mathrm{4}},\end{array}$$

where α is albedo, SW and LW are incoming short- and long-wave radiation, and σ is the Stefan–Boltzmann constant. R_n is partitioned into latent (LE), sensible (H), ground, and canopy heat fluxes. When in equilibrium and over annual and longer time periods, since the ground and canopy do not gain or lose heat systematically, the sum of LE and H fluxes equals net radiation ($R_{\mathrm{n}}=\mathrm{LE}+H$).

Table 1 Plant functional types (PFTs) represented in CTEM and their relation to CLASS PFTs.

CTEM simulates terrestrial ecosystem processes by prognostically tracking carbon in three living vegetation components (leaves, stems, and roots) and two dead carbon pools (litter and soil) for seven non-crop and two crop PFTs that map directly onto the CLASS PFTs (Table 1). The terrestrial ecosystem processes simulated in this study include photosynthesis, autotrophic respiration, heterotrophic respiration, dynamic leaf phenology, and allocation of carbon from leaves to stem and root components. These processes are described in a sequence of papers detailing parameterization of photosynthesis, autotrophic and heterotrophic respiration (Arora, 2003), dynamic root distribution (Arora and Boer, 2003), phenology, carbon allocation, biomass turnover, and conversion of biomass to structural attributes (Arora and Boer, 2005). A full description of CTEM can be found in the appendix of Melton and Arora (2016) for CO₂-related processes and Arora et al. (2018) for CH₄-related processes.

The structure of CTEM is shown in Fig. 1; the original three live vegetation pools (leaves, stem, and roots) are indicated by L, S, and R subscripts, respectively, and the two dead carbon pools (litter or detritus and soil carbon) are indicated by D and H subscripts, respectively). Time varying fluxes in and out of these carbon pools (C_L, C_S, C_R, C_D, and C_H; kg C m⁻²) makes them prognostic variables in the model. The corresponding rate change equations for amount of carbon in the three live vegetation components (leaves, stem, and roots) in the original model version are represented by

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{{\mathrm{dC}}_i}{\mathrm{d}t}}=a_{{\mathrm{f}}_i}\left(G-E_{\mathrm{m}}-E_{\mathrm{g}}\right)-D_i=a_{f_i}N-D_i,\end{array}$$

where the index i corresponds to each of the live vegetation pools ($i=\mathrm{L},\mathrm{S},\mathrm{R}$), $a_{{\mathrm{f}}_i}$ represents allocation fraction for a given vegetation component, E_m is vegetation maintenance respiration, E_g is vegetation growth respiration, and D_i represents the litter loss. G is canopy GPP and the primary photosynthetic carbon flux which drives all other carbon fluxes in and out of the model's various pools. Although photosynthesis is a function of the total leaf area (LA, m²), model parameterizations use LAI (m² m⁻²) since all carbon fluxes and pools are represented per unit of land area.

$N=G-E_{\mathrm{m}}-E_{\mathrm{g}}$ is the canopy NPP and therefore $a_{{\mathrm{f}}_i}N$ represents the fraction of NPP allocated to the three vegetation components. Growth respiration, E_g, is estimated as a fraction of the positive gross canopy photosynthetic rate after maintenance respiration has been accounted for (equation A28; Melton and Arora, 2016). $E_{\mathrm{a}}=E_{\mathrm{m}}+E_{\mathrm{g}}$ is the autotrophic respiration; therefore, $N=G-E_{\mathrm{a}}$. When heterotrophic respiration (E_h) is accounted for, net ecosystem productivity (NEP) is calculated as NEP $=G-E_{\mathrm{a}}-E_{\mathrm{h}}=N-E_{\mathrm{h}}$. Positive NEP values indicate that land is gaining carbon from the atmosphere. Combined autotrophic and heterotrophic respiration (E_a+E_h) are referred to as the ecosystem respiration (E_r).

2.1.1 Addition of NSC pools

For the modifications made in this study, first, NSC pools are included in each of the live vegetation components (leaves, stem, and roots). The total biomass (kg C m⁻²) for each of these components is divided into its non-structural and structural components (indicated by subscripts NS and S) as shown in Fig. 1. ${\mathrm{C}}_{\mathrm{L}}={\mathrm{C}}_{\mathrm{L},\mathrm{NS}}+{\mathrm{C}}_{\mathrm{L},\mathrm{S}}$ and similarly for C_S and C_R. The fraction of NPP allocated to each live vegetation component is first moved to its non-structural part, and a flux of carbon from the non-structural to the structural part provides carbon to the structural part. Once the carbon is moved from non-structural to a structural part of a component it cannot be moved back. Since NPP includes respiratory losses, this essentially implies that respiratory carbon losses are assumed to occur from the non-structural part. This is consistent with observational and other modelling studies (e.g., Hoch et al., 2003; Sperling et al., 2015; Li et al., 2016), which show that plants' NSC stores become depleted due to respiration during the dormant season and environmental stresses associated with cold temperatures and drought conditions. Litter losses in the model, however, occur from both the structural and non-structural parts of leaves and stem and root components.

The modified rate change equations for carbon in the non-structural and structural parts of leaf (Eq. 3), stem, and root (Eq. 4) components are thus written as

$$\begin{array}{}\text{(3)}& {\displaystyle }& {\displaystyle }\begin{array}{l}{\displaystyle \frac{{\mathrm{dC}}_{\mathrm{L},\mathrm{NS}}}{\mathrm{d}t}}=a_{f_{\mathrm{L}}}N-D_{\mathrm{L},\mathrm{NS}}-F_{\mathrm{ns}\mathrm{2}\mathrm{s},\mathrm{L}}+T_{\mathrm{S}}+T_{\mathrm{R}}\\ {\displaystyle \frac{{\mathrm{dC}}_{\mathrm{L},\mathrm{S}}}{\mathrm{d}t}}=F_{\mathrm{ns}\mathrm{2}\mathrm{s},\mathrm{L}}-D_{\mathrm{L},\mathrm{S}}\end{array},\text{(4)}& {\displaystyle }& {\displaystyle }\begin{array}{l}{\displaystyle \frac{{\mathrm{dC}}_{j,\mathrm{NS}}}{\mathrm{d}t}}=a_{{\mathrm{f}}_i}N-D_{j,\mathrm{NS}}-F_{\mathrm{ns}\mathrm{2}\mathrm{s},j}-T_j\\ {\displaystyle \frac{{\mathrm{dC}}_{j,\mathrm{S}}}{\mathrm{d}t}}=F_{\mathrm{ns}\mathrm{2}\mathrm{s},j}-D_{j,\mathrm{S}}\end{array};j=\mathrm{S},\mathrm{R},\end{array}$$

where $F_{\mathrm{ns}\mathrm{2}\mathrm{s},i}\left(i=\mathrm{L},\mathrm{S},\mathrm{R}\right)$ represents carbon flux from the non-structural to structural part of a component (leaf, stem, or root), and $T_j\left(j=\mathrm{S},\mathrm{R}\right)$ represents the reallocation (or transfer) of carbon from stem and root components to leaves during the leaf-out period (see Sect. 2.1.2). Note that there are no autotrophic respiration terms in Eqs. (3) and (4) since they are already included in the term N, the net primary productivity. F_ns2s,i is represented as

$$\begin{array}{}\text{(5)}& F_{\mathrm{ns}\mathrm{2}\mathrm{s},i}={\mathit{\mu }}_i{a}_{{\mathrm{f}}_i}Nmax\left[\mathrm{0},\left(\mathit{\eta }-{\mathit{\eta }}_{i,min}\right)\right],\end{array}$$

where μ_i is a non-dimensional coefficient set to 70. Equation (5) attempts to keep the fraction of non-structural to total carbon in a component ${\mathit{\eta }}_i={\mathrm{C}}_{i,\mathrm{NS}}/{\mathrm{C}}_i$ above its minimum specified value ${\mathit{\eta }}_{i,min}$. During periods of negative NPP, for example, as is the case during winter for cold deciduous trees when they do not have their leaves on, F_ns2s,i is set to zero. This represents the attempts by plants to conserve their NSC pools during a period of no productivity. The amount of carbon in non-structural and structural parts of all vegetation components is shown as time-varying variables and therefore so is the ratio of non-structural to total carbon (η_i). The minimum ratio of non-structural to total carbon in a component (${\mathit{\eta }}_{i,min}$) is specified to be 0.05 for the broadleaf cold deciduous PFT considered here, following Li et al. (2016).

The above modifications made to version 2.1.1 of CTEM in regards to the inclusion of NSC pools allow the movement of non-structural carbohydrates among the model's three live vegetation components, in particular, reallocation of non-structural carbohydrates from stem and root components for leaf out at the onset of spring for the broadleaf cold deciduous tree PFT. Distinction between sugar and starch components of NSC pools is not made. This distinction is not necessary to achieve our primary objective of including NSC pools, which is to address the issue of delayed phenology. In addition to incorporating NSC pools, we also adjust allocation fractions for the leaves and stem and root components after summer solstice in response to day length, and we adjust the lower temperature threshold for leaf litter generation due to cold stress. Deciduousness at high latitudes is determined by both day length and temperature (Xie et al., 2015), and these modifications, discussed below, help to improve simulated leaf phenology.

2.1.2 Reallocation of non-structural carbon during leaf-out period

Leaf phenology in CTEM is represented via four phenological states a plant can be in at any given time (Arora and Boer, 2005). These stages include no leaves or dormant state, maximum leaf growth state, normal growth state, and leaf-fall or harvest state. Depending on their deciduousness, CTEM's nine PFTs (Table 1) may or may not go through these four different leaf phenological states. A broadleaf cold deciduous tree transitions through all four states in a year: leafless/dormant state in winter, maximum growth state (following arrival of favourable climatic conditions in spring when all NPP is allocated to leaves), normal growth state (after reaching a threshold LAI, NPP is dynamically allocated to stem and roots in addition to leaves), and finally the leaf-fall state (triggered by unfavourable environmental conditions and with no carbon allocation to leaves). When all the leaves have been shed, the trees go back into the leafless or dormant state again and the cycle repeats itself in the next year. The phenology module of the CLASS–CTEM model is explained in detail in Melton and Arora (2016).

In the original version of the model, when a plant moves into the maximum leaf growth state all NPP is allocated to leaves until a threshold LAI (L_thrs, m² m⁻²) has been grown. L_thrs is about 40%–50% of the maximum LAI a plant can support depending on its stem and root biomass and based on an allometric relationship between green and woody biomass (Melton and Arora, 2016). In the absence of NSC pools in the original model version, photosynthesis during the early leaf-out period is based on a small imaginary number of leaves (referred to as storage LAI). Once the actual LAI exceeds the storage LAI then photosynthesis is based on the actual LAI. Storage LAI is proportional to a plant's stem and root biomass and was intended as a proxy for the size of NSC pools. However, the rate of photosynthesis from a reasonably apportioned storage LAI is still too small to realistically “push out” leaves at the onset of spring in a time period of about 2 weeks as seen in observations. This limitation in the original model version is overcome here with the inclusion of NSC pools. In the modified version of the model used here, a specified fraction of carbon amount needed to reach the threshold LAI is reallocated ($T_j,j=\mathrm{SR}$) from a plant's stem and root NSC pools to the non-structural part of leaves every day until LAI reaches L_thrs. Note that while this reallocation occurs the leaves are still able to photosynthesize and able to increase their biomass as in the original model version, depending on meteorological conditions. The objective of carbon reallocation from stem and roots to leaves is to accelerate the rate of leaf expansion and LAI increase during leaf onset.

The amount of carbon reallocated (kg C m⁻²) from stem and root components to leaves is given by

$$\begin{array}{}\text{(6)}& {\displaystyle }& {\displaystyle }{T}_j=\mathit{\beta }{\displaystyle \frac{L_{\mathrm{thrs}}}{\mathrm{SLA}}}{f}_j;j=\mathrm{S},\mathrm{R}\text{(7)}& {\displaystyle }& {\displaystyle }{f}_j=\left\{\begin{array}{ll}{\displaystyle \frac{{\mathrm{C}}_{j,\mathrm{NS}}}{{\mathrm{C}}_{\mathrm{S},\mathrm{NS}}+{\mathrm{C}}_{\mathrm{R},\mathrm{NS}}}}& \mathrm{if}{\mathit{\eta }}_j>{\mathit{\eta }}_{j,min}\\ \mathrm{0}& \mathrm{if}{\mathit{\eta }}_j\le {\mathit{\eta }}_{j,min}\end{array}\right\};j=\mathrm{S},\mathrm{R},\end{array}$$

where SLA is the specific leaf area (m² kg⁻¹ C⁻¹), β is the reallocation coefficient set to $\mathrm{6.66}\times {\mathrm{10}}^{-\mathrm{3}}$, and fractions $f_j\left(j=\mathrm{S},\mathrm{R}\right)$ ensure that carbon reallocated from stem and root NSC pools is proportional to the size of their NSC pools. Equation (7) also shows that when the fraction of NSC pool relative to total carbon in a component $\left({\mathit{\eta }}_j={\mathrm{C}}_{j,\mathrm{NS}}/{\mathrm{C}}_j\right)j=\mathrm{SR}$ is equal to or drops below its minimum specific value $\left({\mathit{\eta }}_{j,min}\right)$ then reallocation is stopped. Reallocation is only performed during the leaf-out state when trees are in the maximum leaf growth state.

Figure 2 Latitude dependence factor (Γ) (Eq. 8) for reducing allocation fraction to leaves after summer solstice.

2.1.3 Adjustments to allocation fraction to leaves after the summer solstice

CTEM uses dynamically calculated allocation fractions (Arora and Boer, 2005; Melton and Arora, 2016) for leaves, stem, and roots, which are based on the light, water, and leaf phenological status of vegetation. The allocation to the three live vegetation components is based on assumptions that carbon is preferentially allocated: (1) to roots when soil moisture is limiting, (2) to leaves when LAI is low, and (3) to the stem to increase vegetation height and lateral spread when increasing LAI leads to a decrease in light penetration. These allocation fractions are superseded by three additional rules: (1) all carbon is allocated to leaves at the time of leaf out for cold deciduous tree PFTs to accelerate leaf development, (2) allocation fractions are adjusted when necessary to ensure a tree has enough stem and root biomass to support leaves (to satisfy a structural allometric relationship), and (3) a minimum realistic root-to-shoot ratio is maintained for all PFTs.

The seasonality of globally averaged LAI is dominated by the response of Northern Hemisphere's vegetation to increasing temperatures in spring and decreasing temperatures in autumn. When compared to observation-based estimates of globally averaged LAI, the CLASS–CTEM-simulated LAI starts to increase later in spring, shows a peak in August (compared to July in observation-based estimates), and exhibits a much slower rate of decline after reaching its annual maximum, which typically occurs just after the summer solstice in each hemisphere (Fig. 11 of Anav et al., 2013). This last issue is addressed by reducing the allocation fraction for leaves ($a_{{\mathrm{f}}_{\mathrm{L}}}$) of cold deciduous tree PFTs after summer solstice by multiplication with a day-length-dependent factor (Γ) given by

$$\begin{array}{}\text{(8)}& {\displaystyle }& {\displaystyle }\mathrm{\Gamma }={\left[{\displaystyle \frac{d}{d+\left(d_{max}-d\right)\mathrm{0.5}\left(\tan h\left(\frac{\mathit{\pi }}{\mathrm{180}}\left(\mathrm{20}\mathit{\varphi }-\mathrm{800}\right)\right)+\mathrm{1}\right)}}\right]}^{\mathrm{20}},\text{(9)}& {\displaystyle }& {\displaystyle }d=\mathrm{24}-{\displaystyle \frac{\mathrm{24}}{\mathit{\pi }}}\mathrm{acos}\left[max\left(-\mathrm{1},min\left({\displaystyle \frac{\sin \mathit{\varphi }\sin {\mathit{\delta }}_{\mathrm{c}}}{\cos \mathit{\varphi }\cos {\mathit{\delta }}_{\mathrm{c}}}},\mathrm{1}\right)\right)\right],\end{array}$$

where d is the day length at latitude ϕ (radian), d_max is its maximum value (h), and δ_c (radians) is solar declination. Γ varies between 0 and 1 and its behaviour in Fig. 2 shows how allocation to leaves is reduced at a faster (slower) rate closer to poles (Equator) after summer solstice in the Northern Hemisphere (21 June). Below 30^∘ N in the Northern Hemisphere Eq. (8) yields Γ=1, so allocation fraction for leaves is not modified. Deciduousness due to day length and temperature typically does not occur in tropics where it is primarily controlled by soil moisture. Neither do broadleaf deciduous cold trees typically exist in the tropics. Similar behaviour is obtained for the Southern Hemisphere after 21 December. Since the allocation fractions for leaves and stem and root components should add to 1, a decrease in allocation fraction for leaves implies an increase in allocation fraction for stem and root components in the modified version of the model. The use of summer solstice to initiate changes in plant behaviour is reasonable since summer solstice is a trigger for many plant physiological processes. Adjustments to allocation fraction for leaves after summer solstice have also been made by Gim et al. (2017) (their Eq. 6). Luo et al. (2018) showed that summer solstice marks a seasonal shift in plant growth, leaf physiology, and foliage turnover in temperate and boreal forests.

In the original version of the CLASS–CTEM model, continuous allocation of carbon to leaves up to the time until they are completely shed led to increase in LAI throughout the growing season rather than a near-constant value or slowly decreasing LAI after the summer solstice.

2.1.4 Adjustments to the lower air temperature threshold

The CLASS–CTEM model is able to respond to environmental conditions and to transition between different leaf phenological states (Arora and Boer, 2005). Leaf litter (D_L) generation is caused by normal turnover of leaves as well as drought and cold stresses, which all contribute to LAI seasonality.

$$\begin{array}{}\text{(10)}& D_{\mathrm{L}}=C_{\mathrm{L}}\left[\mathrm{1}-e^{\left(-{\mathrm{\Omega }}_{\mathrm{N}}-{\mathrm{\Omega }}_{\mathrm{C}}-{\mathrm{\Omega }}_{\mathrm{D}}\right)}\right],\end{array}$$

where C_L is the leaf carbon pool and ${\mathrm{\Omega }}_{\mathrm{N},\mathrm{C},\mathrm{D}}$ are the leaf loss rates (day⁻¹) associated with normal (N) turnover of leaves and the cold (C) and drought (D) stresses. The leaf loss rate associated with cold stress (Ω_C) is based on Eqs. (A49)–(A50) of Melton and Arora (2016) (shown below as Eq. 11):

$$\begin{array}{}\text{(11)}& {\mathrm{\Omega }}_{\mathrm{C}}={\mathrm{\Omega }}_{\mathrm{C},\max }{L}_{\mathrm{cold}}^{\mathrm{3}},\end{array}$$

where Ω_C,max is the maximum leaf loss rate due to cold stress and L_cold is a scalar that varies between 0 and 1 as follows

$$\begin{array}{}\text{(12)}& L_{\mathrm{cold}}=\left\{\begin{array}{ll}\mathrm{1}& ,T_{\mathrm{a}}<\left(T_{\mathrm{cold}}^{\mathrm{leaf}}-\mathrm{5}\right)\\ \mathrm{1}-{\displaystyle \frac{T_{\mathrm{a}}-\left(T_{\mathrm{cold}}^{\mathrm{leaf}}-\mathrm{5}\right)}{\mathrm{5}}}& ,\left(T_{\mathrm{cold}}^{\mathrm{leaf}}-\mathrm{5}\right)<T_{\mathrm{a}}<T_{\mathrm{cold}}^{\mathrm{leaf}}\\ \mathrm{0}& ,T_{\mathrm{cold}}^{\mathrm{leaf}}<T_{\mathrm{a}}\end{array}\right\}.\end{array}$$

$T_{\mathrm{cold}}^{\mathrm{leaf}}$ is a PFT-dependent parameter below which a PFT experiences damage to its leaves and this promotes leaf loss due to cold stress in the model.

Figure 3 Location of the three FLUXNET sites chosen in this study to evaluate the changes made to the CLASS–CTEM parameterizations aimed to improve leaf phenology. Figure adapted from Google Maps.

The original version of the model used a $T_{\mathrm{cold}}^{\mathrm{leaf}}$ parameter value of 8 ^∘C throughout the year. In the modified version of the model used here for the broadleaf cold deciduous tree PFT a $T_{\mathrm{cold}}^{\mathrm{leaf}}$ value of 12 ^∘C is used after summer solstice. For the broadleaf cold deciduous tree PFT, leaf out starts in spring, the maximum LAI occurs between July and September (during the Northern Hemisphere's summer), and the leaves are shed between October and November during autumn. Increasing $T_{\mathrm{cold}}^{\mathrm{leaf}}$ leads to more leaf litter generation due to the cold stress in the autumn and moves the descending side of the LAI curve inwards during autumn.

Table 2 The location of FLUXNET sites, primary species that exist at these sites, their soil physical characteristics, mean annual values of primary meteorological variables, and years of data availability.

2.2 Model evaluation and experimental set-up

2.2.1 Description of FLUXNET sites

We evaluate the performance of the original and modified versions of the CLASS–CTEM framework in simulating leaf phenology at three well-studied sites in the eastern United States, which are selected from the FLUXNET network: (1) Harvard Forest (US-Ha1) located at 42.53^∘ N and 72.17^∘ W, (2) Morgan–Monroe State Forest (US-MMS) at 39.32^∘ N and 86.41^∘ W, and (3) the University of Michigan Biological Station (US-UMB) at 45.55^∘ N and 84.71^∘ W. The location of the three FLUXNET sites is shown in Fig. 3. The selected sites meet our requirement of availability of observation-based LAI data (against which our model results can be evaluated) and are primarily characterized by deciduous broadleaf forests although with different species composition. The LAI measurements are based on an LAI-2000 plant canopy analyzer instrument and details are provided in Urbanski et al. (2007), Schmid et al. (2000), and Gough et al. (2008) for the Harvard Forest, Morgan–Monroe, and University of Michigan sites, respectively. The mean annual climate at these sites and their species composition are summarized in Table 2. While these sites differ somewhat in the climate they experience, they share enough commonalities in climate to exhibit similar seasonal dynamics of LAI. Annual precipitation at these temperate locations (US-Ha1, US-MMS, and US-UMB) is 1189, 1083, and 613 mm, with an annual mean temperature of 8.2, 12.4, and 7.2 ^∘C for each site, respectively. These annual averages are based on the half-hourly meteorological data that are used to drive the CLASS–CTEM model for the time period summarized in Table 2.

The US-Ha1 site is owned by Harvard University. Most of its surrounding area was cleared for agriculture during European settlement in 1600–1700. The trees at the site have been regrowing since before 1900 and are currently characterized by predominantly red oak and red maple, with patches of mature hemlock stand and individual white pine. Climate measurements have been made at the Harvard Forest since 1964. The US-MMS site is owned by the Indiana Department of Natural Resources. Many of the trees in the tower footprint are 60–80 years old. Today, the forest is a secondary successional broadleaf forest within the maple–beech to oak–hickory transition zone of the eastern deciduous forest. Finally, the US-UMB site is located within a protected forest owned by the University of Michigan and consists of mid-aged northern hardwoods, conifer understory, aspen, and old-growth hemlock.

The permeable soil depths are specified at 2.5, 2.5, and 2.62 m at the US-Ha1, US-MMS, and US-UMB sites, respectively. Soil texture information was adapted from the global data set of Zobler (1986) and used to specify the percentage of sand and clay in the model's three soil layers as follows. At US-Ha1, the percentages of sand in the first, second, and third soil layers are specified at 68.5 %, 66.5 %, and 72.25 %, and the percentage of clay at 5 %, 5 %, and 4.25 %, respectively. At US-MMS, the percentages of sand in the first, second, and third soil layers are specified at 21 %, 22.5 %, and 30.25 % and the percentage of clay at 21 %, 23 %, and 23.75 %, respectively. At US-UMB, the percentages of sand in the first, second, and third soil layers are specified at 71 %, 72.5 %, and 73.25 % and the percentage of clay at 7 %, 7 %, and 7.75 %, respectively.

2.2.2 CLASS–CTEM simulations

For the three sites investigated here, we have used version 3.6 of CLASS coupled to version 2.1.1 of the CTEM model and made changes mentioned above in Sect. 2.1. Model performance is evaluated for both the modified and original (without NSC pools) versions against available observation-based estimates of LAI and energy and CO₂ fluxes. Simulations were performed for the broadleaf cold deciduous tree PFT with 100 % fractional cover.

Seven meteorological variables are required to drive the CLASS–CTEM model – air temperature, air pressure, wind speed, incoming short-wave radiation, incoming long-wave radiation, precipitation, and specific humidity. FLUXNET's gap-filled meteorological forcing was obtained for each of the three sites. The data were either available at a half-hourly time step or linearly interpolated from hourly to half-hourly resolution. The meteorological data used to drive the model correspond to the period of 1998–2013 for the site in Harvard forest, 1999–2006 for the site in Morgan–Monroe State Forest and 1997–2013 for the site at the University of Michigan Biological Station.

Table 3 Simulated and observation-based turbulent energy and carbon fluxes and LAI at the three FLUXNET sites.

All simulations were forced with meteorological data from their respective FLUXNET sites repeatedly until model carbon pools reached equilibrium and the annually averaged NEP was close to zero. A specified atmospheric CO₂ concentration of 350 ppm is used at all sites for this spin-up. The real-world forests have, of course, experienced a gradual increase in atmospheric CO₂ concentration, changes in climate, and disturbances over their lifetime. Although not perfect, in the absence of full histories of disturbance and meteorological data at these sites this approach still allows comparison of the seasonality of simulated LAI and primary carbon and energy fluxes with observation-based estimates once the model pools reach equilibrium. One caveat is that the modelled vegetation and soil carbon pools cannot be expected to be exactly the same as in the real world but we still expect them to be reasonable. We have chosen to use an atmospheric CO₂ concentration of 350 ppm to spin up the model pools (while the average CO₂ concentration during the first decade of the 21st century was around 380 ppm) because the terrestrial biosphere is not in equilibrium with the atmospheric CO₂ concentration. The disturbance (fire) module was not activated in these simulations. Observation-based LAI measurements were obtained from the AmeriFlux web site (https://ameriflux.lbl.gov, last access: 14 March 2018). Energy and CO₂ fluxes were obtained from the FLUXNET web site (https://fluxnet.fluxdata.org, last access: 14 March 2018).

3 Results

Model performance is evaluated by comparing simulated LAI and CO₂ fluxes of GPP and NEP, which is our primary focus. We also compare radiative energy fluxes of R_n and LE and H with their observation-based estimates from the modified and the original model versions.

Figure 4 Observed and CLASS–CTEM-simulated averaged daily values of (a) LAI (m² m⁻²), (b) GPP (g C m⁻² day⁻¹), (c) NEP (g C m⁻² day⁻¹), and (d) ecosystem respiration (g C m⁻² day⁻¹) for the US-Ha1 (Harvard Forest) FLUXNET site across all years when data are present. Legends also show the mean annual value of the quantity plotted, except for LAI, which is averaged over the growing season. Root-mean-square error (RMSE) and coefficient of determination (R²) are also shown for simulated values when compared to observation-based estimates.

3.1 LAI and land–atmosphere CO₂ fluxes

Figures 4–6 compare simulated values of LAI, GPP, NEP, and E_r from the two model versions with their observation-based estimates at the US-Ha1, US-MMS, and US-UMB FLUXNET sites. Observation-based measurements are shown in black and simulated mean daily values are shown in red (for the original model version indicated as CLASS–CTEM original) and blue (for the modified version, with NSC pools and other changes indicated as CLASS–CTEM modified). Just like simulated values, the observation-based estimates also represent average daily values across all years for which the data were available. The mean annual values of LAI, GPP, E_r, and NEP are also summarized in Table 3. At all sites, when compared to the original version, the modified version of the model shows a phenological shift of about 2 weeks earlier in the year, which is in better agreement with observed LAI transitions (Figs. 4a, 5a, and 6a). The timing of maximum LAI also improves and shows a shift of about 2 months earlier in the year, from late September and early October to late July and early August. The observation-based estimates of LAI suggest the presence of understory vegetation at two of the three FLUXNET sites (the Monroe and the Michigan sites). The CLASS–CTEM modelling framework does not represent any understory vegetation. Despite this, the model still overestimates maximum LAI at all locations and its implications are discussed further down. At all three sites, the inclusion of non-structural carbon pools (Sect. 2.1.1) and other model modifications (Sect. 2.1.2 to 2.1.4) produces a notable improvement in simulated LAI seasonality, especially during canopy development (i.e., spring and early summer) and its autumn decline.

Figure 5 Observed and CLASS–CTEM-simulated averaged daily values of (a) LAI (m² m⁻²), (b) GPP (g C m⁻² day⁻¹), (c) NEP (g C m⁻² day⁻¹), and (d) ecosystem respiration (g C m⁻² day⁻¹) for the US-MMS (Morgan–Monroe State Forest) FLUXNET site across all years when data are present. Legends also show the mean annual value of the quantity plotted, except for LAI, which is averaged over the growing season. Root-mean-square error (RMSE) and coefficient of determination (R²) are also shown for simulated values when compared to observation-based estimates.

The Morgan–Monroe site (Fig. 5) experiences somewhat warmer temperatures than the Harvard and Michigan sites (Figs. 4 and 6) (mean annual temperature at the Morgan–Monroe is about 4 ^∘C higher than at the other two sites; see Table 2) and as a result the growing season is somewhat longer at the Morgan–Monroe site. The model is able to successfully capture this difference amongst the sites. Overall, the simulated GPP and NEP (Figs. 4b, c, 5b, c, 6b, c) compare reasonably well with observations. Improvements in simulated LAI seasonality lead to concomitant improvements in simulated GPP, especially at the ascending side of the plots when the growing season starts. In the original version of the model the increase in GPP at the start of the growing season is delayed due to delayed leaf out. Note that the simulated GPP values compare well with their observation-based estimates despite the higher simulated LAI. Improvements in simulated GPP also lead to improvements in simulated NEP in Figs. 4c, 5c, 6c, and similar to GPP, especially on the ascending side of the plots at the start of the growing season.

Figure 6 Observed and CLASS–CTEM-simulated averaged daily values of (a) LAI (m² m⁻²), (b) GPP (g C m⁻² day⁻¹), (c) NEP (g C m⁻² day⁻¹), and (d) ecosystem respiration (g C m⁻² day⁻¹) for the US-UMB (University of Michigan Biological Reserve) FLUXNET site across all years when data are present. Legends also show the mean annual value of the quantity plotted, except for LAI, which is averaged over the growing season. Root-mean-square error (RMSE) and coefficient of determination (R²) are also shown for simulated values when compared to observation-based estimates.

The comparison with observation-based estimates of LAI and GPP is not completely straightforward since the observation-based estimates of these two quantities are put together by different communities and the fact that GPP is a derived quantity (as opposed to NEP, which is directly observed). As a result, the observation-based estimates of LAI and GPP are not completely consistent with each other. This is seen in Fig. 5 for the US-MMS site where the modified version of the model results in a better match with observation-based estimates of LAI (Fig. 5a), but it shows a bias towards an early increase in GPP (Fig. 5b). In contrast, in Fig. 6, the simulated GPP in the modified version of the model compares better with its observation-based estimate than the LAI. In this respect, NEP provides a better measure to assess model improvement than GPP. In Figs. 4 to 6 improvements in simulated LAI in the modified version of the model are more consistent with improvements in simulated NEP.

The individual contributions of the three modifications, (1) inclusion of NSC pools, (2) reduced allocation of carbon to leaves after summer solstice, and (3) change in parameter value of temperature threshold for leaf litter loss due to cold stress, made to the model for resulting improvements in LAI, GPP, and NEP are shown in the Supplement.

The annual mean of observation-based NEP values (as shown in the figure legends of Figs. 4 to 6) is positive because Northern Hemisphere temperate land is currently a sink of carbon (Myneni et al., 2001). In contrast, the annual mean of simulated NEP values is close to zero by construction because the model was spun up to an equilibrium state. The positive annual mean of observation-based NEP values, compared to simulated NEP values, can manifest in multiple ways – as primarily higher summer values when NEP values are positive (as for the Harvard Forest site), as higher values through the year (as is mostly the case at the Morgan–Monroe site), and as less negative NEP values during the non-growing season when NEP values are negative (as seen at the University of Michigan site). Regardless of this caveat, the inclusion of NSC pools to advance leaf onset and offset times does lead to an improvement in seasonality of simulated NEP values.

Figure 7 CLASS–CTEM (modified version) simulated values of total (a) and non-structural (c) carbohydrate pools (kg C m⁻²). Panel (b) shows the reallocation of carbon from non-structural stem and root pools to leaves during leaf onset in spring and panel (d) shows the carbon flux from non-structural to structural pools for leaf, stem, and root components (g C m⁻² day⁻¹) for the US-Ha1 (Harvard Forest) FLUXNET site. The plots show mean daily values across all years for which the meteorological data were available after the model pools reached equilibrium.

While photosynthesis primarily depends on the current meteorological conditions and LAI amongst other environmental factors (including atmospheric CO₂ concentration), ecosystem respiration (Figs. 4d, 5d, 6d) depends strongly on the vegetation and soil carbon pool sizes. As a result, if simulated vegetation and soil carbon pools are larger or smaller than observation-based estimates, then so would be the respiratory fluxes. Note also that the simulated annual respiratory fluxes are higher than observed at all three sites (Figs. 4d, 5d, 6d, and Table 3). Had the simulated fluxes been lower than what they are now and closer to their observation-based estimates, then the simulated NEP would have been more similar to observations. Nevertheless, the model simulates the seasonality of ecosystem respiratory fluxes reasonably well. In absence of the long-term disturbance history or meteorological data to drive the model with, the current methodology (in which the model is driven repeatedly with the available observed meteorological data) is reasonable and allows us to assess the seasonality of simulated LAI and land–atmosphere CO₂ fluxes – which is the primary objective of our study.

Figure 8 CLASS–CTEM (modified version) simulated values of total (a) and non-structural (c) carbohydrate pools (kg C m⁻²). Panel (b) shows the reallocation of carbon from non-structural stem and root pools to leaves during leaf onset in spring and panel (d) shows the carbon flux from non-structural to structural pools for leaf, stem, and root components (g C m⁻² day⁻¹) for the US-MMS (Morgan–Monroe State Forest) FLUXNET site. The plots show mean daily values across all years for which the meteorological data were available after the model pools reached equilibrium.

3.2 NSC pools

Figures 7–9 evaluate the seasonal cycle of the NSC pools in leaf, stem, and root vegetation components at the US-Ha1, US-MMS, and US-UMB FLUXNET sites, respectively. There are no observation-based estimates of NSC pools available at the three FLUXNET sites. For the broadleaf cold deciduous tree PFT considered here, the stem carbon pool is the largest (and so are its structural and non-structural parts) and the leaf carbon pool is the smallest (Figs. 7a, c, 8a, c, 9a, c). The amount of non-structural carbon reallocated from the stem and root NSC pools to leaves during leaf onset in early spring (see Sect. 2.1.2) is shown in Figs. 7b, 8b, 9b. Figures 7d, 8d, 9d show the seasonality of the carbon flux from the non-structural to the structural part of the leaf, stem, and root components for the three sites. The seasonality of total stem and root carbon pools is driven mostly by the seasonality of their non-structural parts.

Figure 9 CLASS–CTEM (modified version) simulated values of total (a) and non-structural (c) carbohydrate pools (kg C m⁻²). Panel (b) shows the reallocation of carbon from non-structural stem and root pools to leaves during leaf onset in spring and panel (d) shows the carbon flux from non-structural to structural pools for leaf, stem, and root components (g C m⁻² day⁻¹) for the US-UMB (University of Michigan Biological Reserve) FLUXNET site. The plots show mean daily values across all years for which the meteorological data were available after the model pools reached equilibrium.

For the stem and root components, the non-structural parts contribute about 6 %–10 % to the total pool size. During the early leaf-out period when reallocation from stem and root NSC pools to leaves is taking place (Sect. 2.1.2), the stem's NSC pool becomes depleted. This transfer/reallocation stops after a threshold LAI is achieved. The transfer of NSC from stem and root pools to leaves occurs mostly through the stem (see Figs. 7b, 8b, 9b) since its NSC pool is about 3–4 times larger than the root component. The NSC pool for both components reduces during the period when leaves are not present (and GPP is zero) due to respiratory and litter losses. The pools for both stem and root components are replenished later during the growing season when a sufficient number of leaves have been grown and allocation of carbon to stem and root components is restored. This is seen in Figs. 7d, 8d, 9d), which show the flux of carbon from non-structural to structural leaf, stem, and root components. Early on during the growing season, carbon flux from the leaf NSC pool to its structural part is much higher since the model preferably allocates carbon to leaves as discussed in Sect. 2.1.2. After a threshold LAI is reached, carbon is also allocated to stem and root NSC pools, which subsequently start to allocate carbon to their structural pools and the tree biomass continues to increase. At the end of the growing season, when photosynthesis stops, allocation to all three components and the fluxes from NSC to structural parts terminate. During the dormant winter season NSC pools provide for the respiratory costs.

Figure 10 Observed and CLASS–CTEM-simulated daily net radiation (W m⁻²), latent heat flux (W m⁻²), and sensible heat flux (W m⁻²) for the three FLUXNET sites. Legends also show the mean annual value of the quantity plotted. Root-mean-square error (RMSE) and coefficient of determination (R²) are also shown for simulated values when compared to observation-based estimates.

3.3 Energy fluxes

Figure 10 compares observation-based measurements of LE, H, and R_n fluxes at the three FLUXNET sites, with their simulated values from the two model versions. Annual mean values of these observation-based and simulated radiative and turbulent energy fluxes are also summarized in Table 3. Unlike the simulated fluxes, the annual mean sum of the observed LE and H, averaged over the years for which observations are available, is not equal to the observed R_n. This non-closure of the energy budget is seen at all three sites and is a typical characteristic of eddy-covariance-based flux measurements (Gao et al., 2017). The annual energy budget closure is off by 17 % at the University of Michigan Biological Station, by 20 % at the Harvard Forest, and by 30 % at the Morgan–Monroe sites as seen in Table 3. Keeping this caveat in mind, the model captures the seasonality of radiative and turbulent fluxes shown in Fig. 10 reasonably overall, with the exception of late winter and early spring. During this period, as solar radiation increases R_n is underestimated (Fig. 10a–c) until the canopy approaches a full-leaf state and this leads to an underestimation of H (Fig. 10g–i) and overestimation of LE. This may be caused by an overestimation of canopy transmissivity and underestimation of snow and soil masking by leafless forests with increasing solar elevation (recently observed in unpublished simulations with CLASS–CTEM at the Borden forest, Borden, Ontario, Canada) and may also be exacerbated by the lack of representation of a small evergreen needleleaf fraction at US-Ha1 and a conifer understory at US-UMB. LE is apparently overestimated throughout the year at US-MMS but we suspect this reflects a larger underestimation of LE relative to H in the measured fluxes; Oliphant et al. (2004) found that accounting for long sampling tube damping effects on LE resulted in a 16 % improvement in energy balance closure at this site. The change to an earlier leaf phenology in the modified simulations results in a slightly earlier increase in LE in the spring as well as slightly earlier decreases in autumn at US-Ha1 and US-UMB, but differences are much smaller at US-MMS.

In Fig. 10 the changes in LAI, due to modifications made to the original version of the model, do not significantly affect LE fluxes because of two related reasons. First, at mid- to high-latitude locations where soil moisture constraint is not very large, as is the case at the three sites considered in this study, total evapotranspiration (or LE flux) is controlled by available energy. This is the reason for the expected seasonality in LE flux at these sites, which is characterized by higher values during summer and lower values during winter. Second, since the LE flux at these three sites is controlled primarily by available energy, the resulting implication is that if evaporative demand cannot be met by transpiration then it will be met by evaporation from the soil. As a result, changes in LAI do not significantly affect total evapotranspiration but change the partitioning of evapotranspiration flux coming from transpiration, evaporation of intercepted water on canopy leaves, and evaporation from the soil.

4 Discussion and conclusions

The CLASS–CTEM model, similar to other land surface schemes implemented in other Earth system models, is not tuned for any specific location but is expected to behave realistically at all locations. Model processes correspond to generic PFTs, in this case broadleaf cold deciduous trees, and are not meant to represent specific species differences within a PFT. It is nearly impossible, at present, to determine more than 100 parameters for individual species for the CLASS–CTEM model. As a result, while our three chosen sites are characterized by different species (as shown in Table 2), they must be represented by a single set of parameter values which correspond to the broadleaf cold deciduous PFT.

Previous studies using the CLASS–CTEM model in the context of land–atmosphere CO₂ fluxes and simulated carbon pools have evaluated its performance at point (Arora, 2003; Arora and Boer, 2005; Melton et al., 2015), regional (Garnaud et al., 2015; Peng et al., 2014; Arora et al., 2016), and global (Arora and Boer, 2010; Melton and Arora, 2014, 2016) scales. These studies indicate that the model performance is reasonable. CLASS–CTEM also participated in the TRENDY model intercomparison, the result of which contributed to the Global Carbon Project for the years 2016 and 2017 (Le Quéré et al., 2016, 2018). A typical model evaluation exercise at the global and regional scale compares model-simulated geographical distribution of GPP, vegetation biomass, and soil carbon with their respective observation-based estimates. Point-scale studies, however, typically focus on the simulated seasonality of energy and CO₂ fluxes as is the case in this study. Model evaluation exercises not only help in identifying model limitations but also yield opportunities to improve model performance by tuning model parameters.

We chose to perform equilibrium simulations by forcing the model repeatedly with available meteorological data at a specified CO₂ concentration of 350 ppm. In the absence of meteorological data that show a warming trend over the historical period we would not have been able to properly perform a historical simulation. Our past experience shows that steady-state simulations, as opposed to historical transient simulations, allow an easier interpretation of model modifications in the absence of confounding effects of changing climate and increasing CO₂.

Previous evaluations of the CLASS–CTEM model that highlighted its limitation of delayed leaf phenology (e.g., Anav et al., 2013) were the motivation for this study. NSC pools play an important role during leaf onset for broadleaf deciduous cold trees, but also other PFTs, and their effect in the original model was represented using the concept of imaginary leaves whose LAI is assumed to be directly proportional to non-leaf biomass. Here, we have explicitly included NSC pools in the model framework along with some other changes and these modifications do lead to improvement in simulated leaf phenology and concomitant improvements in the simulated seasonal cycle of GPP and NEP. Improvements in simulated energy fluxes are much harder to detect because the observation-based energy fluxes are affected by non-closure of the energy budget but also because LE fluxes are not as strongly dependent on LAI as GPP.

Despite the simulated LAI being higher than observation-based estimates, the simulated GPP, E_r, and NEP compare reasonably with their observation-based estimates. Possible reasons for higher simulated LAI include higher-than-observed allocation to leaf compared to stem and root components and lower-than-observed leaf turnover and/or leaf respiration rates. The model currently uses a maximum carboxylation capacity (V_c,max, i.e., maximum photosynthetic rate) value of 57 µmol CO₂ m⁻² s⁻¹ for broadleaf cold deciduous trees based on Table 3 of Kattge et al. (2009), who derive V_c,max values for major PFTs using more than 700 data estimates. While in the CLASS–CTEM model photosynthesis is also limited by light and transport capacity rates in addition to carboxylation capacity (see Appendix A2 in Melton and Arora, 2016), V_c,max remains a strong parameter and simulated GPP in the model is proportional to V_c,max. While the model-simulated LAI can be lowered by tuning allocation to leaves, leaf turnover, and/or respiration rates specifically for these sites, this would imply using a V_c,max value higher than that suggested by Kattge et al. (2009) to achieve realistic GPP. It is possible that the average V_c,max value derived by Kattge et al. (2009) is not representative of broadleaf cold deciduous trees in the eastern United States. The simulated LAI in the model is the result of multiple model processes interacting with each other. We note this limitation of the model at these locations and plan to address it in the near future. While LAI is an important determinant of model performance, even more important are the land–atmosphere CO₂ fluxes from an Earth system model perspective since it is the CO₂ fluxes which determine the carbon budget of the atmosphere in a fully coupled Earth system model simulation (Arora et al., 2013).

Plants are extremely complex living organisms which respond to the changes in their physical and chemical environmental conditions using a myriad of adaptations. Our limited understanding of these adaptations comes only from empirical observations of their behaviour and measurement of their physical and chemical responses to environmental changes. Models typically represent only a fraction of this understanding because model structures depend on the purpose of the model and the number of details that can be represented reasonably in a model's framework. In hindsight, the omission of NSC pools in the original model version was a structural error and while the conceptual imaginary leaves tried to mimic the fast growth rate of leaves during leaf onset at the arrival of favourable environmental conditions, they were not completely successful in capturing the real-world behaviour. In addition, in a past exercise, we also used higher V_c,max values at the beginning of the growing season to accelerate the rate of growth of leaves (based on Bauerle et al., 2012, and Alton, 2017) but this also did not sufficiently help to address the slower-than-observed rate of growth of LAI at the start of the growing season. Unlike physical models, which describe a physical process, modelling of biological response to changes in environmental conditions is more complex. While there may be underlying physical laws that determine the response of plants to changes in environmental conditions, we can only interpret this from a biological point of view. Dynamic vegetation models and land surface schemes parameterize biological functioning using mathematical formulations to reproduce empirical observations and modellers' conceptual understanding of how the biology works. The inclusion of NSC pools in the CLASS–CTEM framework is based on the same philosophy.

The implementation of NSC pools in the CLASS–CTEM modelling framework presented in this study is meant specifically to address the problem of delayed leaf phenology. NSC pools also play a vital role in the overall health of the plants as mentioned earlier in the introductory section. During periods of limited photosynthesis, trees depend solely on stored NSCs to maintain basic metabolic functions, produce defensive compounds, and retain cell turgor (Sperling et al., 2015). A period of continuous drought, for instance, will gradually reduce the size of NSC pools and this can be used as a trigger to initiate drought-related mortality in the model, or alternatively NSC pools may be used to allow leaf growth during a short-term dry period to represent resilience (Mitchell et al., 2013). The inclusion of NSC pools, also lays the groundwork to implement a nitrogen (N) cycle in the CLASS–CTEM framework since modelling V_c,max as a function of leaf N content requires leaf N content in the non-structural part of the leaves.

In conclusion, modifications to the CLASS–CTEM framework made in this study to address the problem of delayed leaf phenology yield improvements to simulated seasonality of LAI at the three FLUXNET sites considered here. These improvements, especially the inclusion of NSC pools also lay the groundwork for future model development and inclusion of new processes.