A simple time-stepping scheme to simulate leaf area index, phenology, and gross primary production across deciduous broadleaf forests in the eastern United States

Terrestrial plants play a key role in regulating the exchange of energy and materials between the land surface and the atmosphere. Robust models that simulate both leaf dynamics and canopy photosynthesis are required to understand vegetation–climate interactions. This study proposes a simple time-stepping scheme to simulate leaf area index (LAI), phenology, and gross primary production (GPP) when forced with climate variables. The method establishes a linear function between steady-state LAI and the corresponding GPP. The method applies the established function and the MOD17 algorithm to form simultaneous equations, which can be solved together numerically. To account for the time-lagged responses of plant growth to environmental conditions, a time-stepping scheme is developed to simulate the LAI time series based on the solved steady-state LAI. The simulated LAI time series is then used to derive the timing of key phenophases and simulate canopy GPP with the MOD17 algorithm. The developed method is applied to deciduous broadleaf forests in the eastern United States and is found to perform well for simulating canopy LAI and GPP at the site scale as evaluated using both flux tower and satellite data. The method also captures the spatiotemporal variation of vegetation LAI and phenology across the eastern United States compared with satellite observations. The developed time-stepping scheme provides a simplified and improved version of our previous modeling approach to simulate leaf phenology and can potentially be applied at regional to global scales in future studies.

1 Introduction

Terrestrial plants play a key role in regulating the exchange of energy and materials (e.g., radiation, heat and moisture, carbon, and trace gas) between the land surface and the atmosphere (Beer et al., 2010; Zhu et al., 2017). The canopy structure and characteristics govern solar radiation interception and absorption (Ni-Meister et al., 2010; Yuan et al., 2013). Plants control water transpiration and photosynthetic carbon fixation through processes from transient changes in leaf stomatal conductance to seasonal variation in foliage dynamics (Eagleson, 2005). In turn, external environmental conditions, such as sunlight, temperature, and water and nutrient availability, selectively determine plant form and function (Bonan, 2008). Numerical models that integrate multidisciplinary knowledge allow us to understand and predict the interactions between terrestrial ecosystems and the climate.

Developments of terrestrial biosphere models essentially seek accurate solutions to the simulation of energy and material exchanging fluxes between ecosystems and the atmosphere. In terrestrial biosphere models, plant canopies are typically characterized using leaf area index (LAI; leaf area per unit ground area) because plant leaf is the basic organ that intercepts solar radiation for photosynthesis and transpiration (Li et al., 2018; Yuan et al., 2013). The exchanging fluxes of energy and materials over a vegetation canopy can then be modeled as a function of environmental conditions (e.g., sunlight, soil moisture, temperature, and humidity) and vegetation LAI (Ding et al., 2014). The development of satellite remote sensing technology offers large-scale observations for vegetation monitoring and a number of modeling approaches have been developed to quantify and simulate land surface fluxes based on climate variables and satellite-derived LAI. These methods include both light use efficiency models (e.g., the Carnegie–Ames–Stanford approach (CASA) model, Potter et al., 1993; the MOD17 algorithm, Running et al., 2004; the vegetation photosynthesis model (VPM), Xiao et al., 2004; the eddy covariance light use efficiency (EC-LUE) model, Yuan et al., 2010; and the two-leaf light use efficiency (TL-LUE) model, He et al., 2013) and process-based models (e.g., the boreal ecosystem productivity simulator (BEPS; Liu et al., 1997), the Breathing Earth System Simulator (BESS; Ryu et al., 2011), the growing production day model (GPD; Xin, 2016), and the revised Simple Biosphere (SiB2; Sellers et al., 1996b) model). Despite being different from each other in the representation of vegetation processes, these methods have been successfully used for applications from field to global scales. While remotely sensed vegetation data perfectly complement canopy process models, the ability to dynamically simulate vegetation LAI is fundamental to enhancing our abilities to predict terrestrial ecosystem processes.

Modeling vegetation leaf dynamics via climate variables requires in-depth understanding of plant phenological processes. This modeling is still largely empirical to date and contributes considerable uncertainties to current terrestrial biosphere models (Richardson et al., 2012). One common method for simulating vegetation phenology is to predict the timing of key phenophases such as spring onset and autumn senescence in a growing season (Hufkens et al., 2018; Liu et al., 2018). For example, most phenology models originate from the growing degree day (GDD) model, a method first proposed by De Réaumur dating back to 1735 (De Réaumur, 1735). The GDD model assumes that plant leaf onset begins when daily mean temperatures accumulated from a fixed date reach a critical threshold. Studies have identified the fact that various environmental factors other than temperature could affect plant phenology to certain degrees (Polgar and Primack, 2011), and therefore efforts have been made to improve the GDD model by adding different influential factors, such as photoperiod, soil temperature, humidity, and soil moisture (Chuine et al., 1999; Hufkens et al., 2018; Liu et al., 2018; Melaas et al., 2013; Yang et al., 2012). Land surface models like the Community Land Model (Oleson et al., 2013) and the Biome-BGC model (White et al., 2000) use a set of complicated and empirical equations to predict the timing of key phenophases across plant functional types. Arora and Boer (2005) developed a carbon-gain-based scheme that initiates leaf onset when environmental conditions are beneficial for the plant in carbon terms to produce new leaves and initiates leaf offset when environmental conditions are unfavorable with incurred carbon losses for plants. Another method for vegetation phenology modeling is to simulate the entire LAI time series over a growing season. For example, the DeNitrification DeComposition model uses an optimal seasonal growth curve of plant LAI and then calculates environmental stresses of water and nitrogen to limit daily carbon and nitrogen allocation to plant leaves (Yu et al., 2014). The growing season index as proposed by Jolly et al. (2005) is a widely used method that could simulate seasonal phenology curves using photoperiod, air temperature, and vapor pressure deficit. While these studies have greatly benefitted the development of the leaf phenology models, evaluation on 14 land surface models in deciduous forests suggested that almost all models predicted the start of the season earlier or the end of the season later than observations, and the model biases were typically 14 weeks or more. It is therefore necessary to improve the current phenology models.

The physiological processes of leaf phenology and canopy photosynthesis are interrelated. Plants absorb carbon dioxide to accumulate biomass through photosynthesis and then redistribute the photosynthetic gain to organs such as leaves, roots, and stems to optimize carbon gain. Given limited external resources, plants have evolved to effectively allocate photosynthate to organs in response to environmental conditions to maximize photosynthetic carbon gain, the fundamental bioenergy for survival (Givnish, 1986). The strategy of biomass allocation among growth, maintenance, and reproduction in a continuously changing environment directly determines whether plants can persist under natural competition pressure both interspecies and intraspecies (Bonan, 2002). In essence, new leaf phenology models may need to account for the processes of canopy photosynthesis more closely and explicitly than the current leaf phenology models.

Xin (2016) proposed a parameterization scheme to simulate vegetation productivity and phenology simultaneously. The method, named as the growing production day (GPD) model, uses canopy gross primary production (GPP) instead of air temperature as an indicator that synthesizes various environmental factors on plant photosynthesis to track how the environment is suitable for vegetation growth. Analogous to the method that derives reference evapotranspiration, the developed method defines a hypothetic canopy with fixed LAI to model potential GPP under certain environment conditions. Similar to the GDD model, the GPD model predicts vegetation spring onset to occur when the accumulated reference GPP reaches a critical threshold. The method has been successfully applied to the biomes of evergreen needleleaf forest, deciduous broadleaf forest, and grassland. To allow for predicting the entire LAI time series over a growing season, Xin et al. (2018) further improved the GPD model by proposing a linear function between LAI and GPP at the steady state. The proposed function and the canopy GPP model (i.e., modeling GPP as a function of LAI and climate variables) together form a closed system of equations that includes both vegetation GPP and LAI. The improved GPD model uses the numerical approach, a method that gives an initial value and then iterates to the convergence of the solution, to solve the closed system of equations and derive LAI in the steady state. The improved GPD model then applies the simple moving-average method to the steady-state LAI to obtain the modeled LAI time series. The improved method allows for the modeling of LAI time series in addition to the timing of individual phenophases. There remain shortcomings to overcome for broad applications of the GPD model. First, the simple moving-average method, despite being widely used in many studies, is empirical and cannot be used within the framework of models that operate at incremental time steps. Second, the developed GPD model that includes many subtle vegetation processes, such as canopy radiative transfer, leaf stomatal conductance, leaf transpiration, leaf photosynthesis, and soil evaporation, requires various climate input data and is computationally intensive for regional to global applications.

Aiming to solve the abovementioned problems, the objectives of the study are to (1) develop a time-stepping scheme to simulate both leaf dynamics and vegetation productivity and (2) simplify the GPD model to allow for long-term applications at a large scale. Given that the phenology modeling in deciduous broadleaf forest, a biome that has distinct seasonal growing cycles, still has large uncertainties (Melaas et al., 2016), this study chooses to simulate leaf dynamics for the deciduous broadleaf forests across the eastern United States. If successful, such a method can potentially be used for future applications to other biomes.

2 Methods and materials

2.1 Modeling steady-state leaf area index

One difficulty in vegetation phenology modeling is that the timescale associated with leaf allocation far exceeds that of many other vegetation processes. Unlike leaf photosynthesis that approaches equilibrium within 1 min and stomatal functioning that reaches the steady state in minutes (Sellers et al., 1996a), leaf dynamics take days or even months in response to weather variation (Zeng et al., 2013). Xin et al. (2018) first put forward the concept of steady-state leaf area index, i.e., canopy LAI when time approaches infinity while the environmental conditions remain unchanging. An alternative biological explanation to steady-state LAI is the maximum canopy LAI that an environment can sustain infinitely by its own photosynthetic activities. Supposing that the carrying capacity of canopy LAI is proportional to the total canopy photosynthetic rate under a given environment, the steady-state LAI can be modeled as follows:

$$\begin{array}{}\text{(1)}& {\mathrm{LAI}}_{\mathrm{s}}=m{\mathrm{GPP}}_{\mathrm{s}},\end{array}$$

where LAI_s denotes the steady-state leaf area index, and m denotes the constant ratio of steady-state leaf area index to environmental capacity denoted by GPP_s, which is the steady-state gross primary production.

The above equation, despite having a simple form, provides a critical function that complements the canopy photosynthesis model. Only parameter m is dependent on plant functional type and can be quantified from field measurements as the average ratio of LAI to GPP at canopy closure (i.e., the time when both canopy LAI and GPP reach equilibrium). Studies have developed various canopy photosynthesis models, such as light use efficiency models and process-based models. Our previous studies (Xin, 2016; Xin et al., 2018) implemented a sophisticated canopy model that assembles the sub-models of canopy radiative transfer, leaf stomatal conductance, leaf transpiration, soil evaporation, and leaf photosynthesis. Although the method has been successfully applied to different biomes, the model structure is complicated for studies at the regional to global scales. To simulate canopy photosynthesis, this study implements the MOD17 algorithm, a big-leaf light use efficiency model that uses routine satellite products (Running et al., 2004). The use of the MOD17 algorithm could greatly simplify the modeling processes and reduce the required climate variables, thereby allowing for broad applications. A brief description of the MOD17 algorithm is provided here and details can be found from the user guide of the MODIS GPP product (Running and Zhao, 2015).

Based on the MOD17 algorithm, vegetation GPP corresponding to the steady-state leaf area index can be modeled as follows:

$$\begin{array}{}\text{(2)}& {\mathrm{GPP}}_{\mathrm{s}}=\mathrm{PAR}\times {\mathrm{FPAR}}_{\mathrm{s}}\times {\mathit{\epsilon }}_{\max }\times f\left(T\right)\times f\left(\mathrm{VPD}\right),\end{array}$$

where GPP_s denotes the gross primary production corresponding to the steady-state leaf area index; PAR denotes photosynthetically active radiation; FPAR_s denotes the fraction of photosynthetically active radiation corresponding to the steady-state leaf area index; ε_max denotes maximum light use efficiency; and f(T) and f(VPD) denote the scalar functions that account for the limitation of temperature and vapor pressure deficit, respectively, on canopy photosynthesis.

The fraction of photosynthetically active radiation can be modeled as follows (Turner et al., 2006):

$$\begin{array}{}\text{(3)}& {\mathrm{FPAR}}_{\mathrm{s}}=\mathrm{1}-\exp \left(-k{\mathrm{LAI}}_{\mathrm{s}}\right),\end{array}$$

where k denotes the canopy light extinction coefficient and LAI_s denotes the steady-state leaf area index.

The environmental scalars can be modeled as follows.

$$\begin{array}{}\text{(4)}& {\displaystyle }& {\displaystyle }f\left(T\right)=max\left(min\left({\displaystyle \frac{\mathrm{TMIN}-{\mathrm{TMIN}}_{\min }}{{\mathrm{TMIN}}_{\max }-{\mathrm{TMIN}}_{\min }}},\mathrm{1}\right),\mathrm{0}\right)\text{(5)}& {\displaystyle }& {\displaystyle }f\left(\mathrm{VPD}\right)=max\left(min\left(\mathrm{1}-{\displaystyle \frac{\mathrm{VPD}-{\mathrm{VPD}}_{\min }}{{\mathrm{VPD}}_{\max }-{\mathrm{VPD}}_{\min }}},\mathrm{1}\right),\mathrm{0}\right)\end{array}$$

TMIN denotes daily minimum air temperature; TMIN_min and TMIN_max denote the lower and upper thresholds of daily minimum air temperature for vegetation photosynthetic activities, respectively; VPD denotes daily vapor pressure deficit; and VPD_min and VPD_max denote the lower and upper thresholds of daily vapor pressure deficit for vegetation photosynthetic activities, respectively.

Given the environmental conditions (i.e., given daily values of photosynthetically active radiation, minimum air temperature, and vapor pressure deficit), Eqs. (1) and (2) together form simultaneous equations, meaning that there are two unknown variables (i.e., LAI and GPP at the steady state) and two different general equations. One may derive an analytic solution if both equations have simple forms. But because the dependence of GPP on LAI is nonlinear, deriving the analytic solution is complicated and we could apply the numerical approach to obtain the solutions. Because LAI_s increases as a linear function of GPP_s in Eq. (1) and GPP_s increases as a logarithmic function of LAI_s in Eq. (2), the simultaneous equations have one and only one nonzero solution of LAI_s. To obtain the nonzero solution, the numerical approach starts with a guess value of LAI_s and then iterates to obtain the approximated solution of LAI_s until converging. Note that the numerical approach is widely used in land surface models. For example, as the stomatal resistance, CO₂ partial pressure at the leaf surface, internal leaf CO₂ partial pressure, and leaf net photosynthesis are dependent on each other, the Community Land Model 4.5 uses the numerical approach to solve stomatal resistance and leaf photosynthesis iteratively until the internal leaf CO₂ partial pressure converges. For every day, daily photosynthetically active radiation, daily minimum air temperature, and daily vapor pressure deficit are used as forcing data to calculate LAI_s for the corresponding day. Because photosynthetically active radiation, minimum air temperature, and vapor pressure deficit vary throughout the year, the calculated LAI_s vary from day to day.

2.2 Modeling leaf area index, phenology, and gross primary production

Because the physiological processes through which plants allocate photosynthates to leaves do not respond instantaneously to climate variation, there is a need to simulate vegetation LAI as lagging behind the steady state. One method to account for the time-lagging effect is to apply the simple moving-average method to buffer abrupt changes from individual events in the time series. Our previous study applied the simple moving-average method to model LAI as the unweighted mean of the previous LAI_s as follows (Xin et al., 2018):

$$\begin{array}{}\text{(6)}& \mathrm{LAI}(n+\mathrm{1})={\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n{\mathrm{LAI}}_{\mathrm{s}}\left(i\right),\end{array}$$

where LAI(n+1) denotes leaf area index at the n+1 day; n denotes the number of days; i denotes an index starting from 1 to n; and LAI_s denotes the steady-state leaf area index.

The simple moving-average method, while useful in vegetation phenology modeling, is suitable for retrospective analysis rather than prediction, and, importantly, it does not match most land surface models that operate at incremental time steps. Analogous to the method used to simulate leaf stomatal conductance in response to environmental variation, this study proposes a time-stepping scheme to simulate LAI realistically as lagging behind the steady state by a simple restricted growth model (Sellers et al., 1996a) as follows:

$$\begin{array}{}\text{(7)}& {\displaystyle }& {\displaystyle \frac{\mathrm{dLAI}\left(n\right)}{\mathrm{d}t}}=k_{\mathrm{l}}\left[{\mathrm{LAI}}_{\mathrm{s}}\left(n\right)-\mathrm{LAI}\left(n\right)\right],\text{(8)}& {\displaystyle }& {\displaystyle }\mathrm{LAI}(n+\mathrm{1})=\mathrm{LAI}\left(n\right)+{\displaystyle \frac{\mathrm{dLAI}\left(n\right)}{\mathrm{d}t}},\end{array}$$

where t denotes the time; k_l denotes a time constant that reflects the fact that photosynthesis does not instantaneously lead to new or big leaves; and LAI(n) and LAI_s(n) denote the leaf area index and the steady-state leaf area index at the n time step, respectively.

In the time-stepping scheme, vegetation LAI does not change much during winter or summer as the current LAI is close to LAI_s, whereas vegetation LAI increases (or decreases) during spring (or autumn) as the current LAI is less (or greater) than LAI_s. For example, when the environment turns favorable for plant growth in spring, LAI_s exceeds LAI and dLAI∕dt is positive such that the modeled canopy LAI increases. Note that the method developed here essentially uses the canopy photosynthetic capacity (i.e., the steady-state gross primary production) instead of air temperature as a synthesized indicator to track the suitability of the environment for plant growth in a time series, and therefore the developed method is referred to as the simplified growing production day (SGPD) model following our previous studies (Xin et al., 2018).

Given the modeled LAI time series, both vegetation phenology and canopy GPP can be easily modeled (Xin et al., 2018). Various approaches have already been developed to derive the timing of key phenophases, such as spring onset and autumn senescence, from seasonal LAI trajectories. This study models the phenological transition dates using a simple method that derives the first spring and last autumn dates at which LAI reaches 20 %, 50 %, and 80 % of the seasonal amplitudes (Richardson et al., 2012). The selected relative amplitudes (20 %, 50 %, and 80 %) correspond to different plant growth stages over a growing season. Because the MOD17 algorithm only requires LAI, daily minimum temperature, daily vapor pressure deficit, and daily photosynthetically active radiation as model inputs, the canopy GPP is simply modeled by substituting the modeled LAI time series and the climate variables into the MOD17 algorithm. For the first day of spring when the modeled LAI is zero, the modeled fraction of photosynthetically active radiation is zero and the modeled GPP is zero. As times move forward, the modeled LAI increases and the modeled GPP increases but is still dependent on other climate variables such as solar radiation, temperature, and vapor pressure deficit.

2.3 Comparative studies using the growing season index

The growing season index (GSI), a widely used method in vegetation phenology modeling (Jolly et al., 2005), allows for the modeling of seasonal LAI time series rather than individual phenophases and is implemented to make direct comparisons with the SGPD model. The GSI model performs comparably to or even outperforms other terrestrial biosphere models on predicting the timing of key phenophases for deciduous broadleaf forests (Melaas et al., 2013).

The instantaneous GSI is first derived based on the work of Jolly et al. (2005) as follows:

$$\begin{array}{}\text{(9)}& i_{\text{GSI}}=i_{\text{TMIN}}\times i_{\text{VPD}}\times i_{\text{Photo}},\end{array}$$

where i_GSI denotes instantaneous growing season index, and i_TMIN, i_VPD, and i_Photo denote the instantaneous scalar functions that account for the constraints of daily minimum air temperature, vapor pressure deficit, and photoperiod, respectively, on vegetation growth.

The scalar functions for i_TMIN, i_VPD, and i_Photo have mathematic forms similar to Eqs. (4) and (5) and are obtained the same as defined in Jolly et al. (2005) as follows.

$$\begin{array}{}\text{(10)}& {\displaystyle }& {\displaystyle }{i}_{\mathrm{TMIN}}=max\left(min\left({\displaystyle \frac{\mathrm{TMIN}-{\mathrm{TMIN}}_{\min }}{{\mathrm{TMIN}}_{\max }-{\mathrm{TMIN}}_{\min }}},\mathrm{1}\right),\mathrm{0}\right)\text{(11)}& {\displaystyle }& {\displaystyle }{i}_{\mathrm{VPD}}=max\left(min\left(\mathrm{1}-{\displaystyle \frac{\mathrm{VPD}-{\mathrm{VPD}}_{\min }}{{\mathrm{VPD}}_{\max }-{\mathrm{VPD}}_{\min }}},\mathrm{1}\right),\mathrm{0}\right)\text{(12)}& {\displaystyle }& {\displaystyle }{i}_{\mathrm{Photo}}=max\left(min\left({\displaystyle \frac{\mathrm{Photo}-{\mathrm{Photo}}_{\min }}{{\mathrm{Photo}}_{\max }-{\mathrm{Photo}}_{\min }}},\mathrm{1}\right),\mathrm{0}\right)\end{array}$$

TMIN denotes daily minimum temperature; TMIN_min and TMIN_max denote the lower and upper thresholds of daily minimum air temperature for vegetation photosynthetic activities, respectively; VPD denotes daily vapor pressure deficit; VPD_min and VPD_max denote the lower and upper thresholds of daily vapor pressure deficit for vegetation photosynthetic activities, respectively; “Photo” denotes daily photoperiod; and Photo_max and Photo_min denote the lower and upper thresholds of the daily photoperiod for vegetation photosynthetic activities, respectively.

LAI can be modeled as the simple moving average of the instantaneous GSI scaled using maximum LAI as follows:

$$\begin{array}{}\text{(13)}& {\displaystyle }& {\displaystyle }\mathrm{GSI}={\displaystyle \frac{\mathrm{1}}{n_{\mathrm{day}}}}\sum _{i=\mathrm{0}}^{n_{\mathrm{day}}-\mathrm{1}}{i}_{\mathrm{GSI}},\text{(14)}& {\displaystyle }& {\displaystyle }\mathrm{LAI}=\mathrm{GSI}\times {\mathrm{LAI}}_{\max },\end{array}$$

where GSI denotes growing season index at the n day; n_day denotes the number of days; i denotes an index starting from 0 to the previous day; i_GSI denotes the instantaneous growing season index; and LAI_max denotes the maximum leaf area index at canopy closure.

It is noteworthy that the instantaneous GSI uses the product of the scalars of minimum temperature, vapor pressure deficit, and photoperiod as an indicator to track the potential canopy photosynthetic capacities on a daily basis. Both the GSI model and the SGPD model, despite having different forms, share the same modeling idea. To understand the differences between the simple moving-average method and the time-stepping method, the GSI model is also implemented with the simple restricted growth model as follows:

$$\begin{array}{}\text{(15)}& {\displaystyle }& {\displaystyle }{\mathrm{LAI}}_{\mathrm{s}}=i_{\mathrm{GSI}}\times {\mathrm{LAI}}_{\max },\text{(16)}& {\displaystyle }& {\displaystyle \frac{\mathrm{dLAI}}{\mathrm{d}t}}=k_{\mathrm{l}}\left({\mathrm{LAI}}_{\mathrm{s}}-\mathrm{LAI}\right),\end{array}$$

where i_GSI denotes the instantaneous growing season index; LAI_max denotes the maximum leaf area index at canopy closure; k_l denotes a time constant that accounts for the lagged responses of plant leaf allocation to climate variation; and LAI and LAI_s denote the leaf area index and the steady-state leaf area index, respectively.

With the modeled LAI time series, the phenological transition dates are then retrieved based on the seasonal amplitude ratio method, the same way as processing the LAI time series derived from the SGPD model. Vegetation GPP is modeled by substituting the modeled LAI time series into the MOD17 algorithm.

Table 1 Site information for the studied flux towers of deciduous broadleaf forests.

2.4 Model comparison and parameterization

This study compares four different modeling approaches, including the results simulated using both the SGPD model and the simple moving-average method (hereinafter referred to as SGPD-SMA), using both the SGPD model and the time-stepping scheme (hereafter referred to as SGPD-TS), using both the GSI model and the simple moving-average method (hereafter referred to as GSI-SMA), and using both the GSI model and the time-stepping scheme (hereafter referred to as GSI-TS). The commonly used metrics, including the Pearson correlation coefficient (R), the coefficient of determination (R²), the root mean square error (RMSE), and the mean bias error (MBE), are derived for model assessment and comparison.

As the MOD17 algorithm is a well-parameterized model, this study applies model parameters from the literature directly. Following the user guide of the MODIS GPP product (Running and Zhao, 2015), key parameters in the MOD17 algorithm are set as ε_max=1.165 g C MJ⁻¹, ${\mathrm{TMIN}}_{\min }=-\mathrm{6.0}$ ^∘C, TMIN_max=9.94 ^∘C, VPD_min=0.65 kPa, and VPD_max=1.65 kPa. The light extinction coefficient of the canopy is 0.5. The parameter that defines the ratio of leaf area index to environmental capacity is set as m=0.58 m² (leaf area) g C⁻¹ day⁻¹ as quantified using the average ratio of LAI to GPP at canopy closure using flux tower data. The canopy maximum LAI is set as 5.80 based on the maximum 95th percentile of satellite-derived LAI across sites and years (Xin et al., 2018). The parameter n_day in the simple moving-average method and the parameter k_l in the time-stepping method control the response of plant leaf allocation to environmental variation. The parameter n_day is set as 21 days and the parameter k_l is calibrated as 0.080 day⁻¹.

2.5 Study materials and preprocessing

We evaluate our approach at the site scale using both flux tower data and remote sensing data and at the regional scale using both climate data and remote sensing data for deciduous broadleaf forests in the eastern United States. For the site-scale studies, all the flux tower sites of deciduous broadleaf forests (Table 1) that are available on the AmeriFlux website (http://ameriflux.lbl.gov/, last access: 13 January 2019) were used for analysis. As the developed SGPD model is a simplified version of our previous modeling approach, the site-scale modeling studies only require daily incoming solar radiation, minimum air temperature, vapor pressure deficit, photoperiod, LAI, and GPP data. Daily incoming solar radiation, vapor pressure deficit, and GPP have already been provided in the level 4 flux tower data, whereas daily minimum air temperature was processed from the half-hourly gap-filled level 2 data, and daily photoperiod as required by the GSI model was computed based on Eq. (17) as a function of geolocation and the day of the year (Allen et al., 1998). As the MODIS LAI has been found to match field measurements well for deciduous broadleaf forests in the eastern United States (Myneni et al., 2002), the 8-day 500 m MODIS LAI version 6 products (MOD15A2H) that are downloaded from the Land Processes Distributed Active Archive Center (https://lpdaac.usgs.gov/, last access: 13 January 2019) were used as the reference data. Canopy LAI at each site was extracted from MOD15A2H for the pixel that contains the corresponding site. The extracted 8-day MODIS LAI, if identified as poor quality in MOD15A2H, was replaced using the three-point median-value moving-window technique. Spikes in the LAI time series were removed using the Hampel filter and then gap-filled using the autoregressive modeling approach (Akaike, 1969). The obtained 8-day LAI time series were further smoothed using the Savitzky–Golay filter and then linearly interpolated to generate daily time series. The phenological transition dates were retrieved from daily LAI time series using the method that derives the first spring and last autumn dates at which LAI reaches 20 %, 50 %, and 80 % of the seasonal amplitudes, respectively (Richardson et al., 2012).

$$\begin{array}{ll}\text{(17)}& {\displaystyle }& {\displaystyle }\mathrm{Pho}={\displaystyle }& {\displaystyle \frac{\mathrm{24}}{\mathit{\pi }}}\arccos \left(-\tan \left(\mathit{\phi }\right)\tan \left(\mathrm{0.409}\sin \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }}{\mathrm{365}}}\mathrm{DOY}-\mathrm{1.39}\right)\right)\right)\end{array}$$

“Pho” denotes daily photoperiod, φ denotes latitude, and DOY denotes the day of the year.

Our regional-scale studies used both climate data and satellite remote sensing data from 1982 to 2016. The daily 1000 m Daymet version 3 dataset (Thornton et al., 2012) was downloaded from the Oak Ridge National Laboratory Distributed Active Archive Center (http://daymet.ornl.gov/, last access: 13 January 2019). The Daymet dataset provided daily incoming solar radiation, minimum temperature, vapor pressure, and photoperiod data, and we derived daily vapor pressure deficit as the difference between average saturated vapor pressure and vapor pressure. Two different satellite LAI products, including the Global Land Surface Satellite (GLASS) dataset (Xiao et al., 2014) spanning from 1982 to 2014 and the MODIS LAI dataset (Myneni et al., 2002) spanning from 2001 to 2016, were used for the regional studies. The 8-day GLASS LAI product was generated at 0.05^∘ resolution using AVHRR data for the time period from 1982 to 1999 and at 1000 m resolution using MODIS data for the time period from 2000 to 2012. The 8-day satellite LAI data across the eastern United States were processed the same way as the site-scale data to obtain daily LAI time series. Because seasonal LAI amplitudes for each individual pixel could vary from year to year, the 2001–2010 average seasonal LAI amplitudes were used as a baseline to derive the start of the season (SOS) and the end of the season (EOS) for each pixel for each year as the dates when seasonal LAI reaches 50 % of the multi-year average seasonal LAI amplitude. The growing season length (GSL) was derived as the difference between EOS and SOS. A 500 m MODIS-based land cover map was obtained from the USGS Land Cover Institute (https://landcover.usgs.gov/, last access: 1 December 2018). The land cover map was generated by choosing the land cover classification with the highest overall confidence for each pixel in 10-year (2001–2010) Collection 5.1 MODIS land cover type (MCD12Q1) data (Broxton et al., 2014). The 500 m land cover map was resampled to 1000 m resolution using the majority resampling approach and reprojected to a Lambert conformal conic projection to mask areas that are not covered by deciduous broadleaf forests.

3 Results

3.1 Site-scale modeling

Figure 1 shows an example for the simulated time series of LAI and GPP using data acquired at the US-UMB in 2004. The LAI time series simulated using both the SGPD-SMA and SGPD-TS methods are consistent with that obtained from MODIS. The LAI simulated using both the GSI-SMA and GSI-TS methods also captures the observed seasonal variation, but the modeled phenophases obviously have a leading phase in spring and a lagging phase in autumn compared with observations. For both the SGPD model and the GSI model, the results derived using the time-stepping method are consistent with those derived using the simple moving-average method, indicating that the time-stepping method is an effective way to reflect the lagging responses of plant leaf allocation to environmental conditions. By substituting the time series of LAI derived from different modeling approaches into the MOD17 algorithm, all the simulated GPP time series could match the flux tower measurements. Daily fluctuation in the observed GPP time series is largely due to variation in solar radiation from day to day. The GPP modeled using both the GSI-SMA and GSI-TS methods have slight overestimates in the phenological transition periods, like spring and autumn, and match well with the flux tower observations in summer and winter.

Figure 1 The modeled and measured daily time series of (a) leaf area index and (b) gross primary production are shown for the flux tower site of US-UMB in 2004. The reference LAI time series in (a) are derived from the MODIS data and the reference GPP time series in (b) are obtained from flux tower measurements.

Figure 2 shows the regression analysis between the modeled and satellite-derived LAI. Overall, the SGPD model outperforms the GSI model on modeling LAI. When evaluated against the MODIS LAI data, the SGPD-SMA and SGPD-TS models achieved R² of 0.887 and 0.890, respectively, and RMSE of 0.804 and 0.778 m² m⁻², respectively, whereas the GSI-SMA and GSI-TS models achieved R² of 0.746 and 0.759, respectively, and RMSE of 1.356 and 1.303 m² m⁻², respectively. Both the GSI-SMA and GSI-TS models simulate LAI reasonably in summer and winter but overestimate LAI in spring and autumn, and therefore the strong correlations between the GSI-modeled and MODIS-derived LAI are largely due to the underlying seasonality of deciduous broadleaf forests. It is noteworthy that the time-stepping method and the simple moving-average method, despite having different mathematical expressions, generate nearly the same simulation results. The R² values between the SGPD-TS model and the SGPD-SMA model and between the GSI-TS model and the GSI-SMA model are 0.989 and 0.994, respectively, and the regression lines are close to the lines of equity, indicating that the time-stepping method is an alternative representation for the simple moving-average method.

Figure 2 Scatter plots are shown for comparisons (a) between the SGPD-SMA LAI and the MODIS LAI, (b) between the SGPD-TS LAI and the MODIS LAI, (c) between the SGPD-TS LAI and the SGPD-SMA LAI, (d) between the GSI-SMA LAI and the MODIS LAI, (e) between the GSI-TS LAI and the MODIS LAI, and (f) between the GSI-TS LAI and the GSI-SMA LAI on a weekly basis. All available site-year flux tower data were included in the analysis. The solid lines denote the 1 : 1 lines and the dashed lines denote the regression lines.

Table 2 The performance of the modeled timings of phenophases as evaluated against satellite observations. The timings of phenophases were derived based on the dates at which the leaf area index reaches 20 %, 50 %, and 80 % of seasonal amplitude. Positive mean bias error (MBE) indicates that the modeled spring onsets are earlier than the observed ones and negative MBE indicates the opposite.

^* Correlation is significant at the 0.001 level.

Table 2 lists the statistical metrics that illustrate the model performance on predicting the timing of different phenophases. As evaluated against satellite observations, the SGPD-SMA model retrieves the spring onset dates well when LAI reaches 50 % seasonal amplitude and the obtained correlation coefficient is 0.718 with RMSE of 13.04 days. The SGPD-TS model performs comparably to the SGPD-SMA model and the resulting correlation coefficients are all significant except for the dates on which autumn LAI reaches 80 % seasonal amplitudes. The SGPD-based models generally outperform the GSI-based models as the achieved correlation coefficients are higher and the RMSEs are smaller. Both the GSI-SMA and GSI-TS models predict spring onsets earlier than observations by more than 30 days and predict autumn senescence later than observations by more than 20 days. By comparison, the SGPD-TS model predicts the dates on which spring and autumn LAI reaches 50 % seasonal amplitudes well with MBE of only −2.56 and −2.86 days, respectively.

Figure 3 Scatter plots are shown for comparisons (a) between the GPP modeled using SGPD-SMA LAI and the flux tower GPP, (b) between the GPP modeled using SGPD-TS LAI and the flux tower GPP, (c) between the GPP modeled using MODIS LAI and the flux tower GPP, (d) between the GPP modeled using GSI-SMA LAI and the flux tower GPP, and (e) between the GPP modeled using GSI-TS LAI and the flux tower GPP on a weekly basis. All available site-year flux tower data were included in the analysis. All modeled GPP values were derived using the MOD17 algorithm. The solid lines denote the 1 : 1 lines and the dashed lines denote the regression lines.

The modeled and measured GPPs are compared in Fig. 3 to understand the performance of GPP modeling. Compared with the flux tower measurements, the results modeled using SGPD-SMA, SGPD-TS, GSI-SMA, and GSI-TS LAI achieved R² values of 0.768, 0.773, 0.722, and 0.719, respectively, and RMSE values of 2.273, 2.239, 2.577, and 2.535 g C m⁻² day⁻¹, respectively. The modeled results using the GSI-based LAI have higher errors, in terms of both RMSE and MBE, than those using SGPD-based LAI. The accuracies of the modeled GPP using SGPD-based LAI are only slightly lower than using the MODIS-based LAI directly. Because GPP is also a function of a range of climate variables, improvements in modeling LAI do not lead to the same amount of improvements in the modeled GPP. The modeling results obtained based on the simple moving-average method are nearly the same as those obtained based on the time-stepping method. Given the high degrees of consistency between the simple moving-average method and the time-stepping method in modeling LAI, phenology, and GPP, only the results obtained using the time-stepping method are shown and discussed in the regional studies as presented in the following section.

Figure 4 The spatial distributions are shown for (a) the 2001–2010 mean MODIS LAI, (b) the 2001–2010 mean SGPD-TS LAI, (c) the 2001–2010 mean GSI-TS LAI, (d) RMSE between SGPD-TS LAI and MODIS LAI, (e) RMSE between GSI-TS LAI and MODIS LAI, (f) MBE between SGPD-TS LAI and MODIS LAI, and (g) MBE between GSI-TS LAI and MODIS LAI across the eastern United States. The units for both RMSE and MBE are m² (leaf area) per m² (ground area).

3.2 Regional-scale modeling

Figure 4 shows the spatial distributions of the 10-year (2001–2010) mean LAI and associated errors as derived from remote sensing data and model simulations. The SGPD-TS method captured the spatial pattern of the satellite-derived LAI well, including the decreasing gradients from south to north and the decreases in mountain areas (Fig. 4a and b). The 10-year mean LAI derived from the GSI-TS method (Fig. 4c) also shows a decreasing trend from south to north but the modeled LAI is much larger than the MODIS LAI. Because the GSI-TS method defines the maximum leaf area index for the growing season, the overestimation of the modeled 10-year mean LAI is primarily due to model overestimates in the spring and autumn phenological transitions. Compared with MODIS observations, RMSE and MBE obtained by the SGPD-TS method are much smaller and distributed more evenly than those obtained by the GSI-TS method. RMSE for the GSI-TS LAI exhibit a decreasing north–south gradient, implying that the model accuracies are lower in southern areas than in northern areas. MBEs for the GSI-TS model are greater than 0.5 m² m⁻² for most areas. When comparing SGPD-TS LAI with MODIS LAI, RMSEs are less than 0.5 m² m⁻² and MBEs are minor across the study region. The amplitudes of the error metrics in the regional-scale studies are consistent with those in the site-scale studies. Note that some studies applied the multi-year mean LAI as derived from remote sensing data to simulate land surface processes. The results obtained here indicate that the SGPD-TS method can be used alternatively to provide multi-year mean LAI time series via climate variables for land surface studies.

Figure 5 The spatial distributions are shown for (a) the start of the season (SOS) derived from MODIS LAI, (b) SOS derived from SGPD-TS LAI, (c) SOS derived from GSI-TS LAI, (d) the end of the season (EOS) derived from MODIS LAI, (e) EOS derived from SGPD-TS LAI, (f) EOS derived from GSI-TS LAI, (g) the growing season length (GSL) derived from MODIS LAI, (h) GSL derived from SGPD-TS LAI, and (i) GSL derived from GSI-TS LAI using 10-year (2001–2010) mean data across the eastern United States. The embedded subplots show comparisons between modeled and MODIS-derived phenological metrics for SOS, EOS, and GSL, respectively.

The spatial distributions for the 10-year mean phenological metrics, including the start of the season (SOS), the end of the season (EOS), and the growing season length (GSL), are shown in Fig. 5. The SGPD-TS method predicts lower SOS (i.e., earlier spring onset), higher EOS (i.e., later autumn senescence), and longer GSL in southern areas than in northern areas. The spatial distributions of all phenological metrics derived using SGPD-TS LAI agree well with those derived using MODIS LAI. From the statistical analysis as shown in the subplots, the phenological metrics derived from the SGPD-TS method achieved correlation coefficient values of 0.879, 0.552, and 0.844, RMSE values of 8.13, 7.54, and 13.73 days, and MBE values of 0.71, −2.82, and −3.54 days, for SOS, EOS, and GSL, respectively, compared to those derived from MODIS data. Although the spatial distributions of the phenological metrics derived from the GSI-TS method match those derived from the satellite observations, the modeled results have considerable biases, with RMSE values of 38.05, 14.37, and 51.58 days and MBE values of −36.33, 12.91, and 49.23 days for SOS, EOS, and GSL, respectively. Consistent with the site-scale studies, the GSI-TS method predicts spring onset much earlier and autumn senescence later than the satellite-derived data, resulting in a large overestimation of the growing season length.

Figure 6 displays the multi-year phenology anomalies that are spatially averaged for deciduous broadleaf forest across the eastern United States. The use of phenology anomalies relative to the 2001–2010 average instead of absolute values makes the results directly comparable. The SGPD-TS method captured the interannual variation of vegetation phenology retrieved from the remote sensing data. When comparing the SGPD-TS method with MODIS (2001–2016) data, the correlation coefficients are 0.896 (p<0.001), 0.650 (p=0.006), and 0.817 (p<0.001) for SOS, EOS, and GSL, respectively. When comparing the SGPD-TS method with GLASS (1982–2014) data, the correlation coefficients are 0.554 (p=0.001), 0.717 (p<0.001), and 0.637 (p<0.001) for SOS, EOS, and GSL, respectively. The SGPD-TS method outperforms the GSI-TS method on capturing the long-term trends of vegetation phenophases, as the correlation coefficients obtained using the GSI-TS method are lower and sometimes insignificant. Yearly fluctuations in EOS derived using the GSI-TS method are smaller than those derived from both the SGPD-TS method and the satellite data. The correlation coefficients between the GLASS data and the MODIS data for SOS, EOS, and GSL from 2001 to 2014 are 0.892, 0.412, and 0.288, respectively. There are only 14 years of overlap between these two different datasets and the correlations are insignificant for both the derived EOS and GSL. The SOS and EOS derived from the GLASS data have much larger variation in 1982–2000 than in 2001–2010. Note that the 8-day GLASS LAI product was generated at 0.05^∘ resolution using AVHRR data from 1982 to 1999 and at 1000 m resolution using MODIS from 2000 to 2012. The significantly reduced interannual variability for SOS, EOS, and GSL after 2000 in the GLASS data suggests that the use of AVHRR and MODIS data in the GLASS dataset could contribute uncertainties to the satellite-derived phenological metrics. Both Figs. 5 and 6 indicate that the SGPD-TS method is reliable in capturing the spatiotemporal patterns of regional vegetation phenophases.

Figure 6 The spatially averaged phenology anomalies relative to the 2001–2010 average are shown for (a) the start of the season (SOS), (b) the end of the season (EOS), and (c) the growing season length (GSL). SOS and EOS are derived as the date on which LAI first and last reaches 50 % of the seasonal amplitudes and GSL is derived as the difference between EOS and SOS. The shaded areas denote the standard deviation of the corresponding phenophases across spaces.

Figure 7 compares the simulated GPP using the MOD17 algorithm and LAI derived from different approaches. The 10-year average annual GPP obtained using SGPD-TS LAI has a similar spatial pattern to that obtained using MODIS LAI and has lower values than that obtained using GSI-TS LAI. Taking the GPP simulated using MODIS LAI as a reference, the results simulated using SGPD-TS LAI achieve a correlation coefficient of 0.898 with RMSE of 78.78 g C m⁻² year⁻¹ and MBE of 12.22 g C m⁻² year⁻¹, whereas the results simulated using GSI-TS LAI achieve a correlation coefficient of 0.898 with RMSE of 173.45 g C m⁻² year⁻¹ and MBE of 153.43 g C m⁻² year⁻¹. Although the obtained correlation coefficients are close, the SGPD-TS method results in regression lines closer to the 1 : 1 line with smaller bias errors than the GSI-TS method. The zonally average profiles of the 2001–2010 average annual GPP as shown in Fig. 7d suggest that the results obtained from the SGPD-TS method are close to those obtained using MODIS LAI, whereas the results obtained from the GSI-TS method have positive biases of approximately 120–180 g C m⁻² year⁻¹ (roughly 10–15 %) across latitudes. Note that the MOD17 algorithm has a positive MBE of 0.247 g C m⁻² day⁻¹ and 0.571 g C m⁻² day⁻¹ when using SGPD-TS LAI and GSI-TS LAI, respectively, as model input data in the site-scale study. The differences in MBE between the two modeling methods are 0.324 g C m⁻² day⁻¹ (or 118.26 g C m⁻² year⁻¹ in equivalence) for the site-scale studies, which are consistent with the regional-scale studies.

Figure 7 Comparisons are shown for (a) the spatial distribution of annual GPP modeled using MODIS LAI, (b) the spatial distribution of annual GPP modeled using SGPD-TS LAI, (c) the spatial distribution of annual GPP modeled using GSI-TS LAI, (d) the zonally averaged profiles of annual gross primary production modeled using LAI derived from different approaches, (e) the regression between GPP modeled using SGPD-TS LAI and MODIS LAI, and (f) the regression between GPP modeled using GSI-TS LAI and MODIS LAI. The simulated daily GPP were first summed for each individual year and were then averaged across years to derive the 2001–2010 average annual GPP as shown in (a–c). The shaded areas in (d) mark the range of the standard deviation. All pixels of deciduous broadleaf forest across the eastern United States are included in the analysis in (e) and (f).

4 Discussion

Here we provide a simple time-stepping solution that allows for the simulation of canopy photosynthesis, leaf area index, and leaf phenology simultaneously. The developed method first proposes a linear function between canopy photosynthetic capacity and steady-state LAI to complement the canopy photosynthesis model and then applies a simple restricted growth model to account for the lagged responses of plant leaf allocation to the natural environment. In essence, the developed method, although having a simple form, has synthesized the impacts of various climate factors on leaf dynamics because any climate variable that influences vegetation photosynthesis would affect the process of plant leaf allocation in the models as well. Consistent with field observations, the simulated LAI increases as the environmental conditions turn favorable for photosynthetic activities such as increases in photoperiod and temperature.

Figure 8 further illustrates the relationship between mean LAI and different variables on a monthly basis. All data were averaged to the monthly timescale such that canopy LAI can be considered as nearly steady state. On a monthly basis, mean LAI has a strong near-linear relationship with mean GPP (R²=0.888) and the slope for the regression without intercept is 0.580, the same as we used in the model simulation. On a monthly basis, mean LAI is strongly correlated with mean temperature (R²=0.799), indicating that temperature is the dominate factor that determines vegetation phenology. Factors like vapor pressure deficit and photoperiod also have positive relationships with mean LAI on a monthly basis. Figure 8 suggests that LAI has a stronger correlation with GPP than with temperature on a monthly basis. Our modeling approach that links canopy GPP with LAI reflects the empirically positive relationship found in Figure 8a.

Figure 8 Scatter plots are shown for the relationship (a) between mean leaf area index and mean gross primary production, (b) between mean leaf area index and mean temperature, (c) between mean leaf area index and mean vapor pressure deficit, and (d) mean leaf area index and mean photoperiod on a monthly basis. All available site-year flux tower data were included in the analysis. All data were averaged to the monthly timescale for analysis, making the point numbers different from the analysis at the weekly timescale in Fig. 2. The dashed lines denote the regression lines. Panel (a) uses the regression without intercept.

The performance of our developed method is largely dependent on the canopy photosynthesis model used. In our previous studies, we developed a process-based canopy photosynthesis model that synthesizes sub-models such as canopy radiative transfer, leaf transpiration, leaf stomatal conductance, leaf photosynthesis, and soil evaporation and applied it for modeling LAI time series. When applying the simple moving-average method, implementing the process-based model in Xin et al. (2018) achieved higher accuracies than implementing the MOD17 algorithm in modeling canopy GPP and LAI as reflected by higher R² and lower errors. The MOD17 algorithm only assumes a monotonic relationship between air temperature and photosynthesis and between vapor pressure deficit and photosynthesis. It also does not account for the impacts of CO₂ on photosynthesis. The use of the MOD17 algorithm in this study thus has limitations in the model structure. It implies that LAI modeling in our developed method will likely benefit from improvements to the canopy photosynthesis model. This study chooses the MOD17 algorithm instead of the sophisticated process-based model because the MOD17 algorithm is well parameterized across biomes and requires quite limited model inputs of climate variables. Successful implementation with the MOD17 algorithm allows for the extension of the developed method to applications across biomes at regional to global scales.

Land surface models that predict vegetation GPP require either satellite-derived LAI input data or the phenology sub-model. The main idea for this study is to improve phenology modeling by providing time series of LAI simulated using climate variables, hence enabling the simulation of GPP forced only by climate variables. Because we implement the MOD17 algorithm instead of the sophisticated process-based model for the purpose of simplicity, one should not expect the GPP simulated based on model-simulated LAI to be more accurate than GPP simulated based on satellite-derived LAI.

The time-stepping scheme developed here is also an improvement over the simple moving-average method as used in our previous studies. The results obtained using the time-stepping method are consistent with the simple moving-average method at the site scale and are shown to be reasonable at the regional scale. Compared to the simple moving-average method, the time-stepping method could be used in models that operate at incremental time steps. For land surface models that include canopy photosynthesis sub-models, the developed method can be embedded into these models as an alternative phenology model if replacing the MOD17 approach with the canopy photosynthesis sub-model. Compared to simple light use efficiency models like the MOD17 algorithm, implementation of the developed time-stepping scheme in land surface models relies on supercomputing for global applications. To better understand the performance of the developed method, one study is now being undertaken to implement the developed method with the Common Land Model for simulating multi-decadal LAI and GPP for global biomes forced only by climate variables.

Applying the developed method to other biomes and other regions still has issues to be solved appropriately. The time-stepping method uses the parameter k_l to account for the time lags of leaf allocation in response to environmental changes. For deciduous broadleaf forests, a biome with strong seasonality, the developed scheme achieved reasonable results with appropriate parameterization. Short vegetation like grasslands tends to respond quickly to abrupt environment changes like precipitation, and tropical ecosystems have strong resilience to short-term environmental variation (Levine et al., 2016; Shen et al., 2011). Another issue is to find the appropriate values of m for different biomes. One way to determine the values of m is to find the regression slope between leaf area index and gross primary production on a monthly basis. Model parameterization, however, requires broad tests. These understandings from observational studies imply that biomes have varied response speeds to the environment, and proper model calibration and assessment are required for the developed method. Using observation data from remote sensing alone is inadequate for model development as satellite-derived LAI could have large uncertainties for some specific biomes other than deciduous broadleaf forests. Fortunately, global flux tower networks and regional phenology observation networks are now established and offer abundant data for comprehensive model assessment.

5 Conclusions

Numerical models provide a basic tool for understanding the interactions between the land surface and the atmosphere. To provide a complete solution to the simulation of plant leaf dynamics and canopy photosynthesis, this study establishes a linear relationship between the steady-state leaf area index and the corresponding canopy photosynthetic capacity. The proposed leaf allocation function complements the canopy photosynthesis model of the MOD17 algorithm to form simultaneous equations that can be solved using the numerical approach. To account for the time lagging of plant leaf allocation in response to climate variation, a time-stepping scheme based on a simple restricted growth model is applied to the solved steady-state leaf area index to obtain time series of leaf area index. The developed method could perform reasonably well on simulating leaf area index, phenology, and gross primary production for deciduous broadleaf forests across the eastern United States over years, as found in both the site-scale and regional-scale modeling studies. Compared to the simple moving-average method, the time-stepping scheme developed here is consistent with and can potentially be embedded into models that operate at incremental time steps. The developed method allows for the simulation of leaf area index and gross primary production simultaneously and provides a simplified and improved version of our previous model as a basis for global applications in future studies.