Implementation of dynamic crop growth processes into a land surface model : evaluation of energy , water and carbon fluxes under corn and soybean rotation

Worldwide expansion of agriculture is impacting the earth’s climate by altering carbon, water, and energy fluxes, but the climate in turn is impacting crop production. To study this two-way interaction and its impact on seasonal dynamics of carbon, water, and energy fluxes, we implemented dynamic crop growth processes into a land surface model, the Integrated Science Assessment Model (ISAM). In particular, we implemented crop-specific phenology schemes and dynamic carbon allocation schemes. These schemes account for light, water, and nutrient stresses while allocating the assimilated carbon to leaf, root, stem, and grain pools. The dynamic vegetation structure simulation better captured the seasonal variability in leaf area index (LAI), canopy height, and root depth. We further implemented dynamic root distribution processes in soil layers, which better simulated the root response of soil water uptake and transpiration. Observational data for LAI, aboveand belowground biomass, and carbon, water, and energy fluxes were compiled from two AmeriFlux sites, Mead, NE, and Bondville, IL, USA, to calibrate and evaluate the model performance. For the purposes of calibration and evaluation, we use a corn–soybean (C4–C3) rotation system over the period 2001–2004. The calibrated model was able to capture the diurnal and seasonal patterns of carbon assimilation and water and energy fluxes for the corn–soybean rotation system at these two sites. Specifically, the calculated gross primary production (GPP), net radiation fluxes at the top of the canopy, and latent heat fluxes compared well with observations. The largest bias in model results was in sensible heat flux (SH) for corn and soybean at both sites. The dynamic crop growth simulation better captured the seasonal variability in carbon and energy fluxes relative to the static simulation implemented in the original version of ISAM. Especially, with dynamic carbon allocation and root distribution processes, the model’s simulated GPP and latent heat flux (LH) were in much better agreement with observational data than for the static root distribution simulation. Modeled latent heat based on dynamic growth processes increased by 12–27% during the growing season at both sites, leading to an improvement in modeled GPP by 13–61 % compared to the estimates based on the original version of the ISAM.


Introduction
Increasing global food demand accelerates deforestation in areas suitable for modern agriculture.Today, croplands and pastures have become the two largest terrestrial biomes, accounting for about 40 % of the planet's land surface (Foley et al., 2005).Additionally, demand for biofuels might accelerate the expansion of croplands in the coming decades.In 2004, about 1 % of global cropland was being used for biofuels, and this share might increase 3-4-fold by 2030 (FAO, 2008).
This rapid transformation of landscape can impact the climate by altering carbon, water, and energy fluxes (Sellers, 1992;McGuire et al., 2001;Matthews et al., 2004;Sitch et al., 2005;Brovkin et al., 2006;Bonan, 2008).While the climate is affected by the expansion of agriculture land, climate change simultaneously affects agriculture.Many crops show Y. Song et al.: Implementation of dynamic crop growth processes into a land surface model positive responses to elevated carbon dioxide and low levels of warming, but higher levels of warming often negatively affect growth and yield (Hatfield et al., 2008;Kucharik and Serbin, 2008;Urban et al., 2012) The overall aim of this study is to evaluate dynamic crop growth processes in a land surface model, the Integrated Science Assessment Model (ISAM) (Jain et al., 2009;Yang et al., 2009;El-Masri et al., 2013), to understand and address the interactions between C3/C4 crop growth and seasonal dynamics of carbon, water, and energy fluxes.Our implementation focuses on the corn-soybean rotation, which is the most common crop rotation practice around the world (Nafziger, 2012).
A number of land surface models incorporate advanced representations of croplands to simulate the relationship better between crop production, land surface characteristics, and energy and water cycles (Tsvetsinskaya et al., 2001;Kucharik and Brye, 2003;Gervois et al., 2004;Bondeau et al., 2007;Osborne et al., 2007;Lokupitiya et al., 2009;Van den Hoof et al., 2011).Tsvetsinskaya et al. (2001) made the first attempt to integrate a corn simulation model into a physical and soil hydrological model, Biosphere-Atmosphere Transfer Scheme -BATS (Dickinson et al., 1993).The coupled model was able to capture the seasonal change in leaf area index (LAI) for corn, and the results demonstrate its importance for the calculation of the surface fluxes of heat, moisture, and momentum.The IBIS was extended to include crops (Donner and Kucharik, 2003) and validated and applied to simulate crop yields, water and energy balance and impacts of agricultural management (Kucharik, 2003;Sacks and Kucharik, 2011).Gervois et al. (2004) implemented a crop simulation model (STICS) (Brission et al., 2002) in the ORCHIDEE land surface model (Krinner et al., 2005) to simulate winter wheat and maize specifically at two sites in western Europe and two sites in the US.Bondeau et al. (2007) have implemented a dynamic representation of carbon allocation, phenology and management practices for a number of crops into Lund-Potsdam-Jena (LPJ) Dynamic Global Vegetation Models (DGVMs) (Sitch et al., 2003).Lokupitiya et al. (2009) developed crop-specific phenology and coupled them to SiB (Sellers et al., 1996a, b).Van den Hoof et al. (2011) implemented dynamic crop growth structure and phenology into JULES-SUCROS (Cox et al., 1999) to study the impact of interactive effect of wheat structure and phenology on land-atmosphere interactions.Most recently, the carbon allocation and phenology algorithms for corn, soybean, and temperate cereals of the Agro-IBIS model (Kucharik and Byre, 2003) have been introduced into community land models (CLMs) (Lawrence et al., 2012) to examine the effects of managed crops on the climate (Levis et al., 2012).
This paper builds upon and extends the approaches of the studies discussed above into ISAM, which has been extensively used in various model inter-comparison studies (Hunzinger et al., 2012;Richardson et al., 2012;Schaeffer et al., 2012;De Goncalves et al., 2013;Kauwe et al., 2013).
While we use a similar carbon assimilation, energy, and hydrological modeling approach, we implement new algorithms to simulate the following processes: (i) crop growth and biomass allocation in five phenology stages, distributing assimilated carbon among above-and belowground parts depending upon both accumulated heat and resource availability, such as light, water, and nutrients (e.g., nitrogen); (ii) development of vegetation structure (LAI, canopy height, and root depth) calculation based on accumulated carbon mass in leaf, stem, and root pools; (iii) vertical and horizontal root growth in soil layers in response to available soil moisture; and (iv) different abscission rates for fresh and old dead leaves.
The newly implemented algorithms, which are described in detail in Sect.2, differ in many ways from the algorithms considered in previous crop growth modeling studies, some of which are discussed above.For example, while the classification of different phenology stages and LAI in the ISAM and most of the other models discussed above are determined according to the fraction of accumulated heat units, the extended version of the ISAM also considers additional phenology stages, such as silk emergence.This accounts for the impact of water stress on crop yield at a critical stage for grain production.Agro-IBIS, SiB-CROP, and the standard ISAM estimate LAI based on leaf carbon by multiplying leaf biomass carbon by specific leaf area (SLA).The extended ISAM further makes a distinction between the green LAI and the standing dead LAI in LAI simulation.The standard ISAM and SiB-Crop simulate the variation of the allocation fractions with cumulated heat units.Adding to this, the extended ISAM further simulates the responses of the allocation fractions to other environmental factors, including water, light, and nutrient availability.While the extended ISAM adopts the STICS-ORCHIDEE and JULES-SUCROS algorithms to calculate the root growth and canopy height, it also calculates vertical and horizontal root growth in soil layers in response to available soil moisture.Overall, unlike crop simulation schemes in other land surface models discussed above, the dynamic crop growth processes implemented in the extended ISAM account for the coupling between carbon biomass dynamics of leaf, stem, root, and grain and vegetation structure (LAI, canopy height and root depth and distribution), as well as environmental factors' (temperature, water, light, nutrients) variability.
Following the implementation of the new processes, the model parameters were calibrated, and model performance was evaluated using observational data (LAI, biomass, and carbon, water, and energy fluxes) from two AmeriFlux sites (Mead, NE, and Bondville, IL) under a corn-soybean rotation system.

Model description
ISAM is a coupled biogeochemical and biogeophysical model with 0.5 • × 0.5 • spatial resolution and multi-temporal resolution from 30 min to 1 yr (Jain et al., 2009;Yang et al., 2009;El-Masri et al., 2013).Each grid cell is occupied by a combination of fractional vegetation, bare soil, and glacier (Meiyappan and Jain, 2012).Here we add two crop functional types (corn and soybean) into the model.The model is driven by the following climate variables at hourly/halfhourly time step: mean surface air temperature, precipitation rate, incoming shortwave radiation, long-wave radiation, wind speed, and specific humidity.There are sunlit and shaded canopies, 10 hydrological and thermal active soil layers, 5 hydrologically inactive and thermally active bedrock layers, 7 vegetation pools, and 8 litter and soil organic matter (SOM) pools in the ISAM (Yang et al., 2009;El-Masri et al., 2013).Carbon assimilation and heat and water fluxes are calculated through coupled canopy photosynthesis and energy and hydrological processes.Carbon assimilation is allocated into vegetation, litter, and soil organic matter (SOM) pools.The C cycle is then coupled with the complete N cycle.The N cycle model accounts for major N processes, including N deposition, N fixation, N mineralization, N immobilization, nitrification, denitrification, and leaching (Yang et al., 2009).
The model variables, parameters, and equations are given in Tables A1 and A2 in Appendix.

Coupled canopy photosynthesis, energy and hydrological balance processes
Carbon assimilation rates and energy and water fluxes are calculated by coupling a leaf temperature, photosynthesis, and stomatal conductance model (Dai et al., 2004;Chen et al., 2010) with an energy and hydrological balance model (Dai et al., 2004;Oleson et al., 2004Oleson et al., , 2008)).The carbon assimilation model is composed of the C3 photosynthesis model (Farquhar et al., 1980;Collatz et al., 1991), the C4 photosynthesis model (Collatz et al., 1992), and the Ball-Berry stomatal conductance model (Ball et al., 1987;Collatz et al., 1991).The stomatal conductance is calculated as a function of net carbon assimilation rate, relative humidity, and CO 2 concentration at the leaf surface.The C3 carbon assimilation rate is co-limited by light availability, RuBisCO efficiency, and carbon compound export ability.The C4 carbon assimilation rate is co-limited by light availability, RuBisCO efficiency, and PEP-carboxylase availability.
The CO 2 compensation point for the C3 biome in the original ISAM is calculated as a function of O 2 partial pressure and temperature-dependent RuBisCO specificity for CO 2 relative to O 2 (Dai et al., 2003).However, this method underestimates the compensation point during the beginning and end of the growing season, resulting in higher gross primary production (GPP) than observed during these two stages of the growing season.Studies (Smith et al., 1976;Kennedy and Johnson, 1981) suggest that the compensation point for young leaves is higher and decreases when the leaves grow, stays constant after their maturity, and increases again during senescence.Following Smith et al. (1976), we calculate the rate of change of the compensation point as a function of leaf age.
Temperature regulates carbon assimilation processes by multiplying temperature functions (Dai et al., 2003) with the maximum carboxylation rate at the reference temperature of 25 • C(V cmax25 ).The effect of soil water availability on carbon assimilation is dependent on V cmax25 and dark respiration and minimum stomatal conductance (Oleson et al. 2008).Moreover, seasonal variation in V cmax25 is calculated based on a day length factor (Bonan et al., 2011).
Leaf level photosynthesis and stomatal conductance are separately scaled to the canopy level for sunlit and shaded leaves by using sun/shade canopy LAI fractions and scaling parameters to represent extinction of nitrogen and light through the vertical canopy (Dai et al., 2004).Carbon cycle equations in the current ISAM are documented in detail in El-Masri et al. (2013).
Energy conservation of the soil-vegetation system in ISAM is calculated as the balance of absorption of net shortwave and long-wave radiation (Rn) by sunlit/shaded canopy and the ground, and emissions of sensible (SH) and latent heat (LH) fluxes from leaves and ground and soil heat fluxes (G).The net solar radiation is calculated by two-stream approximation (Sellers et al., 1996a), which dynamically calculates the interception, reflectance, transmission, and absorption of direct and diffuse radiation by sunlit/shaded canopy and the soil (Dai et al., 2004).The treatment of diffuse radiation in the "two-stream" scheme is based on Bonan et al. (2011) in order to reduce biases in shaded leaf photosynthesis.Vegetation optical characteristics (leaf/stem reflectivity and transmissivity, Table A1), canopy structure (expressed as leaf angle distribution, Table A1) and density (expressed as leaf area index (LAI) and stem area index (SAI)) dynamically control partitioning of canopy-intercepted radiation and ground-intercepted radiation, as well as partitioning of vegetation adsorbed net radiation between sunlit and shaded canopy.Sunlit canopy intercepts direct and diffuse radiation, whereas shaded canopy intercepts only diffuse radiation.
Latent heat transfer to the atmosphere is resolved using canopy transpiration, canopy evaporation from the intercepted precipitation water, condensation of evaporated water, dew formation, and ground evaporation; sensible heat is partitioned into ground and canopy components (Dai et al., 2003).ISAM also considers an additional soil resistance (Sellers et al., 1992) and litter resistance (Sakaguchi et al., 2009) to the humidity transfer from ground to atmosphere.Convergence of aerodynamic properties from thick/thin canopies to that of the ground is ensured based on Zeng and Wang (2007).Surface albedo is resolved into ground albedo (a function of soil color and wetness), exposed vegetation albedo (a function of leaf orientation, leaf/stem reflectivity and transmissivity, and ground albedo), and snow albedo (Zeng et al., 2002).
The hydrological cycle is coupled with the energy cycle through latent heat of evaporation from wet canopy and the ground, transpiration from dry canopy, and watercontent-adjusted soil heat conductivity.Canopy interception and throughfall of precipitation, infiltration, redistribution of soil water within the soil column, surface runoff, and subsurface percolation are calculated using the formulations of Oleson et al. (2004Oleson et al. ( , 2008)).Canopy interception of precipitation and dew formation, as a function of canopy density, determines the dry/wet ratio of the canopy and thus partitioning of evaporation and transpiration from leaves.These processes are coupled with a dynamic root distribution algorithm (Arora and Boer, 2003) and determine the soil water availability for root uptake and transpiration.The soil heat is modeled based on Fick's equation (Dai et al., 2003), whereas the soil water flux is implemented based on Richard's equations (Oleson et al., 2004).The thermal and hydrological properties for each soil layer are estimated based on soil liquid and ice water contents, soil temperature, soil texture and soil organic carbon (SOC), and gravel content (Lawrence et al., 2008).

Implementation of dynamic crop growth processes in the ISAM
The ISAM, as described by El-Masri et al. (2013), is extended to enable the explicit study of dynamic crop growth processes, specifically accounting for the effects of light, water, and nutrient stresses on C3 and C4 crop growth and water and energy fluxes under the soybean-corn rotation system.In particular, we implement crop-specific phenology schemes, dynamic carbon allocation processes, and dynamic vegetation structure growth processes (LAI, canopy height and root depth and distribution).In the following section, we describe each individual dynamic process as it has been implemented in the ISAM for the current study.

Phenology development
The crop phenology begins with the planting of seeds and ends with grain harvest.In between, the phenology is divided into five growth stages: emergence period, initial vegetative period, normal vegetative period, initial reproductive period, and post-reproductive period.The planting date is determined when the following three conditions are satisfied simultaneously (Eq.A1): (1) the mean daily air temperature of the past seven consecutive days is greater than the base temperature (T base ); (2) the mean daily soil temperature of the past seven consecutive days is greater than the crop-specific critical soil temperature for emergence (T soil critical ); and (3) the accumulated growing de-gree days above 0 • C is greater than the crop-specific minimum value (GDD0 min ) (Eq.A1).At the time of planting, seeding rate is given as an input parameter based on field crop management.After planting, the transition of the different growth stages of phenology is determined by the heat unit index (HUI) and the accumulated days for each growth stage (Eqs.A2-A7).The HUI is 0 at planting time and 1 when the crop matures.The required heat value and the total numbers of days for each growth stage are attained from published studies (Darby and Lauer, 2000;McWilliams et al., 2004;USDA-NASS, 2009;USDA-OCE, 2010).These values are further calibrated based on multiyear LAI from the Mead, NE, AmeriFlux site (Verma et al., 2005).Moreover, the cropspecific maximum LAI (LAI max ) is used as a threshold value to control crop growth development (Kim and Wang, 2005).When modeled LAI becomes larger than the threshold LAI, the modeled phenology is transitioned from vegetative to reproductive growth stage, and the leaves begin to turn brown and start falling (Eqs.A3-A7).The model also accounts for extreme cold and warm temperatures on crop yields.The effect of extreme cold temperature on yield, which is referred to here as the frost damage condition, is accounted for by assuming 100 % loss of yield.This condition is activated any time after emergence stage when the mean daily temperature for five consecutive days falls below 273.2K (Darby and Lauer, 2000) (Eq.A8).
As a plant with separate male and female flowering parts, the ear represents the female flower of the corn plant.The silks are the functional stigmas of a corn plant, which collect pollen and transmit the male genetic material to ova and produce viable kernels.Silk emergence from the ear shoot is a critical process in the production of corn grain (Aldrich et al., 1986).When severe drought stress occurs, silk emergence is delayed during the reproductive period.This effect severely decreases corn yields.This effect is triggered in the model when the following two conditions are met simultaneously (Eq.A9): (1) the mean daily temperature of three consecutive days exceeds 303.2K (Shaw, 1988;Rattalino Edreira and Otegui, 2012), and (2) the mean water stress index of three consecutive days is lower than 0.5.Finally, the crops are harvested when they mature (i.e., HUI = 1.0).

Carbon allocation
Assimilated carbon in leaves is allocated to stems, roots, and grain.The leaf component is divided into photosynthetically active (green leaves) and dead (senescent) leaves.Initial carbon is determined based on the amount of carbon stored in the seeds (Eq.A10).During the emergence time, carbon stored in the seed is allocated to the leaf and root based on thermal conditions (Eq.A11).The carbon assimilation through photosynthesis allocates carbon to each vegetation pool (leaf, stem, root, and grain) (Eqs.A18-21).Part of the assimilated carbon is lost through respiration (Eq.A12).Maintenance respiration for each vege-tation pool is calculated as a function of carbon amount, the C : N ratio, and temperature-dependent respiration coefficients (Eq.A13) (Sitch et al., 2003).Temperature-dependent respiration coefficients are calculated based on a specified respiration rate at 20 • C and a Q 10 temperature function (Arora, 2003) (Eq.A14).The Q 10 values for leaf, stem, and root respiration are calculated as a function of leaf, stem, and soil maintenance temperature, respectively (Arora et al., 2003) (Eq.A15).The growth respiration is assumed to be 25 % of the remainder after removing maintenance respiration from GPP (Eq.A16).The partitioning of the growth respiration into each vegetation pool follows the fraction of carbon in each vegetation pool (Eq.A17).
The net assimilated carbon (GPP minus maintenance and growth respirations) allocated to leaf, stem, root, and grain pools is a dynamic process based on temperature, water availability, light, and N to alter the carbon allocation fractions dynamically at each model time step (Eq.A22).The objective of this allocation scheme is that the allocation of carbon into leaves, stems, and roots is adjusted to minimize adverse effects of limited availability of light, water, and mineral nutrients.Accordingly, more carbon is allocated to roots when soil moisture and mineral N are limiting, or it is allocated to stem and leaf when more leaves results in a decrease in light penetrating the canopy (Arora and Boer, 2005;Salter et al., 2003).This dynamic allocation approach is similar to that of Friedlingstein et al. (1999) and Arora and Boer (2005), except that carbon allocation factors for the different stages of phenology vary with HUI (Eqs.A18-26) based on Penning de Vries et al. (1989).
During the vegetative period, allocated carbon in the green leaf and stem increases with HUI in order for the canopy to build and capture increasing amounts of radiation (Eqs.A19 and A25).No new carbon is allocated to the green leaf pool in the initial and post-reproductive period.Instead, the leaf pool loses carbon through maintenance respiration (Eq.A26), and the carbon is allocated to grain with increasing HUI to increase grain filling during the initial reproductive period (Eq.A26).However, no C is allocated to corn grain if the silk emergence is delayed under drought conditions.
The transition from vegetative to reproductive period initiates the process of leaf senescence.During this period, green leaf carbon is reduced, leading to reduced photosynthetic carbon fixation.The conversion of green leaf to dead leaf carbon occurs at death rates that vary due to drought or cold conditions following the formulations of Arora and Boer (2005) (Eqs.A27-30).
During the post-reproductive period, assimilated carbon is only allocated to grain and root pools (Eq.A21).If no green leaves exist before crop maturity, carbon stored in roots and stems is partly reallocated to the grain pool to enhance the grain filling.In order to account for the effect of water stress on grain filling, the reallocation fraction factor is downscaled (Eq.A26).
Finally, a dynamic allocation factor for each vegetation pool is modified to satisfy two conditions (Arora and Boer, 2005).The first is that there must be enough root and stem biomass to support leaf biomass (Eq.A31) The second condition is that a minimum root/shoot ratio must be available to maintain the structure of each crop type (Eq.A32).If the first condition is not satisfied, the carbon is allocated to root and stem.If the second condition is not satisfied, carbon is allocated to root.A fraction of the carbon allocated to the vegetation pools can be lost as litter.Following Arora and Boer (2005), conversion of the root and stem carbon to litter occurs at a fixed turnover rate (Eqs.A33-34).Conversion of dead leaves to litter occurs as a function of fresh dead leaves and accumulated dead leaves produced in previous time steps (Eq.A35).

Calculation of LAI, canopy height and root depth
Total LAI in the model is the sum of green and standing dead LAI, which is calculated as a function of total leaf carbon and specific leaf area (SLA) (Eq.A36).The green LAI is calculated by multiplying the green leaf carbon masses by SLA (Eq.A37), and the standing dead LAI is calculated by subtracting green LAI from the total LAI (Eq.A38).Canopy height, which is used to parameterize atmospheric turbulence above the canopy in the model, is calculated by scaling the maximum canopy height (H max ) with the accumulated aboveground biomass (Arora and Boer, 2005) (Eq.A39).Canopy height increases from 0 to H max with increased aboveground biomass.Root depth and root distribution in each soil layer vary temporally and spatially with accumulation of root biomass (Arora and Boer, 2003) (Eqs.A40-43).The parameter α appearing in Eqs.(A41-42), which ranges from 0 to 1, determines the rate at which root density varies horizontally, and the root depth grows vertically with increased root biomass (Arora and Boer, 2003).As α approaches 1, the more the roots tend to grow vertically.The parameter D norm_profile appearing in Eqs.(A41-42) determines the root distribution under no water stress.Since the allocation of assimilated carbon to root is sensitive to soil water availability (Eq.A21), ISAM-simulated root growth and distribution in each soil layer are dynamically sensitive to soil water availability within each soil layer according to Eqs. (A41-42).The reduced soil water content in the root zone induces water stress, which leads to increased carbon allocation to roots and thus rapid increasing of root biomass according to Eq. (A21).Following Eqs.(A41-42), both root depth and root density in each soil layer increase with increased root biomass.

Description of the site data
Field data for corn and soybeans from two AmeriFlux eddy covariance flux tower sites, Mead, NE, and Bondville, IL, are used to evaluate the performance of the ISAM.Both sites have similar annual mean temperature averaged for the time period 2001-2004, which is around 284 K at Mead, NE, and 285 K at Bondville, IL.However, the annually accumulated precipitation averaged for the time period 2001-2004 at the Bondville, IL, site (about 787 mm) is about 183 mm higher than that at the Mead, NE, site.About 60 mm of this difference is observed during June and July, when precipitation is positively correlated with both corn and soybean yields.In addition, there are differences in soil characteristics at two sites.The Mead, NE, site sits on deep silt clay loam (Suyker et al., 2004), whereas the Bondville, IL, site sits on silt loam (Hollinger et al., 2005).The soils at the Mead site have lower water infiltration rates and lower plant-available water storage abilities than the soil at the Bondville site.The Mead and Bondville sites have been planted with corn and soybean in rotation since 2001 and 1996, respectively.Corn is grown in odd years and soybeans are grown in even years.Weeds are controlled with herbicides, but no tillage or irrigation is used at either site (Meyers and Hollinger, 2004;Suyker and Verma, 2009).The hourly measured carbon, heat, and water exchanges between the atmosphere and canopy, and biweekly measured LAI, leaf carbon, biomass, and annual yield (ftp://cdiac.ornl.gov/pub/ameriflux/data/Level2/Sites_ByName/Mead_Rainfed/) at the Mead rainfed site, Nebraska (41.18 • N, 96.44 • W) (Suyker et al., 2004), are used to calibrate the processes and parameters of the extended version of the ISAM.Then we use the calibrated parameters along with AmeriFlux data from Bondville, Illinois (40.00 • N, 88.29 • W) (Hollinger et al., 2005) (ftp://cdiac.ornl.gov/pub/ameriflux/data/Level2/Sites_ByName/Bondville/) to evaluate the model's performance for carbon (GPP) and energy fluxes (net radiation (Rn) at the top of the canopy, latent heat (LH) and sensible heat (SH) fluxes) between the atmosphere and canopy at both the diurnal and seasonal scale, and seasonal LAI.

ISAM calibration and evaluation
Calibrations of dynamic crop phenology, carbon allocation, and vegetation structure growth (LAI, canopy height and root depth and distribution) processes are performed in five steps.First, the daily LAI data are calculated by interpolating biweekly measured LAI data.Second, the model is run with prescribed daily LAI to calibrate the initial carbon fraction allocated to the leaf, stem, root, and grain (Table 1).This is achieved by comparing observed and measured leaf car-bon biomass, aboveground biomass, and grain yield.Third, instead of using prescribed observed LAI, the model simulates the daily LAI.We then calibrate the parameters that are used in the dynamic phenology simulations (Table 1) by comparing modeled LAI and measured LAI data.These parameters are especially important for capturing seasonal variability in LAI and thus carbon and energy exchange between the canopy and atmosphere.Fourth, the canopy height equation parameter, m (Table 1), is calibrated by comparing simulated and measured canopy height.Finally, the parameters used for the calculations of the dynamic root growth and distribution (Table 1), which have the strongest effect on both carbon and energy fluxes simulations under water stress condition, are calibrated by comparing the observed and calculated root biomass distribution.Since there is not much information available in literature about the root biomass distributions for corn and soybean for the growing seasons studied here, here we use corn root profiles measured for three specific dates in 1980 at the Mead site (Newell and Wilhelm, 1987) to calibrate the root growth direction parameter (α) and root distribution parameter (D norm_profile ) for corn.Due to the lack of site-specific climate forcing data in 1980, we use 1980 NLDAS-2 climate forcing data (Mitchell et al., 2004) to drive the model for this calibration.All other information, such as management seeding rate, planting time, etc., is taken from Newell and Wilhelm (1987).For soybeans, we calibrate α and D norm_profile by comparing measured and modeled soil water content.Calibration is performed by minimizing the total sum of the squares of the difference between simulated and observed data for corn and soybean at the Mead, NE, site.This is realized through automatic optimization using PEST, which is a nonlinear parameter optimization program and can be used with any model (Doherty, 2005).

Model experiments
The site-specific climate and soil data are used in the following model experiments.The climate data for each site are obtained from the AmeriFlux database.The soil texture data for each site are attained from Web Soil Survey (http://websoilsurvey.nrcs.usda.gov/).We spin up the model for each site with a corn-soybean rotation under repeating site climate data from 2001 to 2004 for about 200 yr until the soil temperature and moisture reach a steady state.Then, we run the model with site-specific planting and harvest times from 2001 to 2004 to calibrate and evaluate the model performance.Due to a lack of measured energy balance closure at many sites (Wilson et al., 2002), we perform the energy balance closure correction according to Twine et al. (2000), which preserves the Bowen ratio: where f is the correction factor.N is the total number of available data points at hourly time intervals over two grow- 1.0, 1.0 1.0, 9.5 k2 v2 1.0, 1.0 2.4, 0.0 ing seasons for each crop.Thus, f shows an overall evaluation of energy balance closure over two growing seasons for each crop.The corrected LH or SH is calculated by multiplying the measured LH and SH fluxes with f .All the energy flux terms, except for the storage energy term (S), are measured at the two sites.We assume S for the Bondville site to be 14 and 8 % of hourly Rn during the morning time (7:00-12:00 UTC) of growing seasons for corn and soybean, respectively (Meyers and Hollinger, 2004).This fraction of hourly S gradually reduces to 2 and 0 % of hourly Rn by 17:00 UTC.For the Mead site, Suyker and Verma (2010) have estimated the corrected energy fluxes for the period 2001-2006, which we apply here.

Statistical analysis
The continued hourly/half-hourly observed fluxes have nonrandom errors and biases (Williams et al., 2009).Therefore, regression analysis is not an optimal way to analyze the model performance.Instead, we use the refined Willmott's index (Willmott et al., 2011) method to quantify the degree to which observed hourly GPP and energy and water fluxes are captured by the model.The range of refined Willmott's index, dr, is from −1 to 1.A dr of 1 indicates perfect agree-ment between model and observation, and a dr of −1 indicates either lack of agreement between the model and observation or insufficient variation in observations to test the model adequately.The refined Willmott's index is calculated as Here P i and O i are the individual modeled and observed data, respectively.Ō is the mean of the observed values.N is the number of the paired observation and modeled data.The 2) represents the sum of modeled error magnitude, and the part (Eq.2) represents the sum of the perfect modeled deviation and observed deviation (Willmott et al., 2011).The Willmott index is a more advanced method to evaluate the land surface model performance than previously reported methods (e.g., Medlyn et al., 2005).Some of the statistical methods widely used to evaluate model performance with observed data are Pearson correlation coefficient (r), coefficient of determination (r 2 ), mean absolute error (MAE), root mean square error (RMSE), and others.These traditional methods, however, are not always optimal for evaluating the model-data agreement or disagreement.For example, r or r 2 methods can indicate the overall linear covariation between data and model results, but need to combine with the slope and intercept of the linear regression to evaluate the degree to which the observed results is captured by the model.Contrastingly, Willmott's index is sensitive to differences between the measured and modeled values and itself can express the degree to how much measured variation can be captured by the model (Willmott, 1981).MAE and RMSE are dimensional measures of disagreement, and thus are not independent of data scale and unit.However, dr is a standardized measure of the model disagreement.It is able to calculate the difference between the magnitude of the mean model bias and the observed deviation.The Willmott's index is similar to the model efficiency (ME), which can also estimate the proportion of model bias to measured deviation.However, dr is more natural measure of mean model bias than ME.Unlike ME, which expresses the model bias as the sum of squared differences between modeled and observed data and thus may upscale the modeled biases, dr expresses the model bias as the sum of absolute value of differences between modeled and observed data (Willmott and Matsuura, 2005).Another advantage of the refined Willmott index is that it is bounded on both the upper and lower ends.The refined index has an easily interpretable lower limit of −1.0 and an upper limit of 1.0, so the range of the index is doubled (Willmott et al., 2011).Many other existing indices, including the original Willmott index (Willmott, 1981), are bounded at the top (usually by 1.0) but sometime lack a finite lower bound, which makes assessments and comparisons of poorly performing models difficult.
Here we calculated dr for hourly observed and modeled data, dr h , to examine the degree to which the model represents the hourly variation in observed values.For the Bondville site, the half-hourly observed and modeled data were synthesized into hourly data and then used to calculate dr h .We also calculated dr for daily mean observed and modeled data, dr d , to examine the model performance at a daily timescale.The comparison of dr h and dr d allowed us to evaluate model biases at these different timescales.
In addition, instantaneous soil moisture measurements at the AmeriFlux sites were used to evaluate the modeled soil moisture.

Best fit model results for the calibration site
Figure 1a-f show the model-calibrated results for LAI, aboveground biomass, and canopy height over the period 2001-2004 at the Mead, NE, site.These figures suggest that the calibrated model is able to simulate dynamic phenology development, carbon allocation, LAI and canopy height growth processes over multiyear growing seasons at the Mead, NE, site.The model also captures well the measured trends of root growth with soil depth during the growing season of 1980 for corn at the Mead site (ISAM-Dynamic case results in Fig. 4a-c).Modeled best fit results with measurements for soil water content during the growing seasons for both corn and soybean years are shown in Fig. 3e and j.
Table 2 shows the statistical analysis, and Fig. 2a-h and Fig. 3a-j compare modeled and measured data for GPP and energy fluxes (Rn, LH and SH) over the period [2001][2002][2003][2004].The statistical analysis and direct model-data comparison results suggest that model's estimated carbon assimilation and energy and water fluxes, with the exception of sensible heat flux, are in good agreement with observations at the Mead site.The relatively low dr h and dr d values are found for SH under corn and soybean rotation, suggesting that modeled results are not consistent with observations.The possible reasons for the differences between the modeled and measured data at the Mead site are discussed together with the differences at the Bondville site in the Sect.4.1.3.

ISAM results at the evaluation site
Overall, the model's estimated results for LAI, aboveground biomass, root biomass, and canopy height over the time period 2001-2004 (Fig. 1j, k and l) compare well with corresponding measured values at the Bondville, IL, site, with a few exceptions.The model slightly overestimates aboveground biomass (Fig. 1k) and canopy height (Fig. 1l) for soybeans during the 2004 normal vegetative phenology stage (Julian days from 170 to 200).These results indicate that the model-calibrated parameters are not only able to simulate the dynamic LAI, phenology, carbon allocation, and canopy height processes at Mead, but also able to capture seasonal variability in LAI, biomass and canopy height growth at the Bondville site.The model is also able to capture the basic daily variation in soil moisture for corn and soybean rotations during the growing period 2001-2004 at the Bondville site (Fig. 3o, t), suggesting that the model parameterization for corn and soybean root profiles, which is calibrated based on the Mead site field data, can also be applied to other sites or regions.As it does at the Mead site, the model captures well both measured diurnal and seasonal variability in GPP, net radiation, and latent heat, but not sensible heat flux (Table 2, Figs.2i-p and 3k-t).The possible reasons for these differences are discussed below.

Biases in the model's estimated carbon and energy fluxes GPP
The dr h for GPP under corn and soybean rotation varies from 0.82 to 0.86 at both sites (Table 2), indicating that modelestimated hourly GPP variations for most cases are consistent with the observations (Fig. 2a, e, i, m).The model's results for soybean GPP are improved by regulating the CO 2 compensation point with the leaf age.To examine this effect, we perform an extra experiment (Model − LACO 2 ), which does not consider the effect of leaf age on the CO 2 compensation point, and compare its results with the experiment with the consideration of the effect (Model + LACO 2 ).The comparison (Fig. 4) shows that the implementation of the leaf age effect on the CO 2 compensation point effectively reduces simulated soybean GPP not only at the calibrated Mead site but also at the Bondville site.These downscaled values are in much better agreement with the measured values during the leaf expansion period.The dr d values for GPP under corn and soybean rotation vary between 0.71 and 0.92, suggesting that the model-estimated daily GPP for most of the cases is consistent with the observations (Fig. 3a

SH and LH Fluxes
The dr h values for SH fluxes are 0.68 for both corn and soybean at the Mead site and 0.60 and 0.69 at the Bondville site, whereas the dr h values for LH for corn-soybean at the Mead and Bondville sites vary between 0.84-0.86 and 0.83-0.84,respectively.These results suggest that the model is able to capture most of the variations in observed hourly LH, but has apparent model biases in hourly SH at both sites.The model overestimates SH during the morning hours (UTC 6:00-10:00 a.m.), but slightly underestimates SH during the afternoon hours (after UTC 2:00 p.m.) (Fig. 2c, g, k and o) at both sites.Similar errors in modeled LH are observed at the Mead site (Fig. 2d and h).These discrepancies result from smaller model biases in Rn (Fig. 2b, f, j, and n).The overestimated Rn speeds up the penetration of the stable stratified canopy atmosphere during the morning hours and warms the canopy quickly, leading to a sudden increase in SH and LH fluxes after sunrise.The biases in modeled SH are also observed during the night hours at the Bondville site when the model usually simulates negative SH, instead of the mean zero value of SH in the measurement.The negative modeled SH indicates the simulation of stable stratified atmospheric layers during the nighttime.It is also important to note that observed fluxes through the eddy covariance technique are    usually unreliable during night hours (Goulden et al., 1996), which adds to discrepancies between modeled and observed values.
Compared to modeled diurnal patterns in energy fluxes, the model is better able to capture seasonal patterns in energy fluxes (Fig. 3c-d and r-s), as indicated by the higher dr d values for SH and LH fluxes than dr h values for corn at the Mead site and for soybean at the Bondville site (Table 2).However, model biases in SH and LH are observed during a specific time period for soybeans at Mead and for corn at Bondville, as indicated by the dr d values being lower than the dr h values at the two sites (Table 2).The overestimated SH and underestimated LH are observed during the normal vegetative period (between Julian day 192 and 223) and the initial reproductive period (between Julian day 224 and 252) in 2002 at the Mead site (Fig. 3h).This model discrepancy results from underestimated soil water content during dry periods (Fig. 3j), which reduces the water availability for evapotranspiration, leading to an underestimation of LH flux and overestimation of SH flux.A similar partitioning discrepancy between SH and LH is observed during the normal vegetative period (between Julian day 170 and 200) of the 2001 corn growing season and at the end of the 2003 corn growing season in Bondville (Fig. 3m-n).In addition, the overestimated SH is also partly attributed to the mismatch in energy partitioning between the soil and atmosphere.We find that the model underestimates ground heat flux, leading to an overestimation of SH fluxes (not shown).
Besides uncertainty in modeled fluxes, uncertainties in measurement, such as uncertainties from the measuring equipment, source heterogeneity, and the turbulent nature of the transport process (Richardson et al., 2006), can also contribute to the discrepancy between simulated and mea-sured energy fluxes.Richardson et al. (2006) estimate that overall random measurement error at the Mead, NE, site averaged about 15.5 W m −2 for SH and LH fluxes for the 2002 and 2003 growing seasons.The estimated root mean squared error (RMSE) for the modeled SH and LH fluxes for the Mead site averaged 20.3 and 16.7 W m −2 for 2002 and 2003 growing seasons, which is slightly higher than the measurement uncertainty.This result suggests that current estimates of overall model biases in partitioning SH and LH may be slightly overestimated without considering the measurement uncertainty.In addition, as pointed out by Wilson et al. (2002), the Bowen ratio method (Eq. 1) might have overlooked the biases in the half-hourly/hourly data, such as the tendency to overestimate positive fluxes during daytime and underestimate negative fluxes during nighttime.
The simulated results for carbon and energy fluxes for corn and soybeans at both sites suggest that the model-calibrated parameters can not only be used to simulate corn and soybean growth at water stressed sites accurately, like the Mead site, but can also accurately simulate corn and soybean growth at sites in the normal non-irrigation region, like the Bondville site.Since the extended version of the ISAM couples the dynamic carbon allocation processes with the vegetation structure simulation (LAI, root depth, and distribution at each soil layer), the model has the advantage of simulating the crop yield and resultant carbon and energy fluxes under different environmental conditions and variability.However, the ability of the model to simulate corn and soybean growth on a large scale still needs to be evaluated and also compared with simulations from other land surface models in future studies.

The effects of different dynamic processes on modeled results
In this section we evaluate the importance of four dynamic processes considered in this study -( 1) dynamic carbon allocation, (2) dynamic LAI, (3) dynamic root distribution, and (4) dynamic scale height -by performing the following additional model simulations (Table 3).ISAM-Static: this model is based on fixed carbon allocation, prescribed LAI, prescribed canopy height, as well as prescribed root depth and root allocation factions in each soil layer.All these four processes have been included in the original version of ISAM (El-Masri et al., 2013).
ISAM-StaticC: this is the same as the ISAM-Dynamic experiment, but the carbon allocation parameterization is based on fixed carbon allocation scheme as assumed in the original version of the ISAM.
ISAM-StaticLAI: this is the same as the ISAM-Dynamic experiment, but uses prescribed LAI development as assumed in the original version of the ISAM.
ISAM-StaticR: this is the same as the ISAM-Dynamic experiment, but uses pre-determined root depth and root fraction for each soil layer in space and time as assumed in the original version of the ISAM.
ISAM-StaticH: this is the same as the ISAM-Dynamic experiment, but uses fixed canopy height parameterization as assumed in the original version of the ISAM.
In the original version of the ISAM (El-Masri et al., 2013), referred to here as ISAM-Static, the carbon allocation fractions for leaf, stem, root, and grain pools for each phenology stage are assumed to be the same values as in the case of ISAM-Dynamic but without accounting for limitation of water, light, and nutrients (Table A1), and these fraction values are assumed to be the same for each model year run.The LAI is not dependent on the carbon allocation simulation as in the case of the ISAM-Dynamic experiment; rather the LAI values in the original version of ISAM are attained from multiyear average site-specific MODIS land product subsets (ORNL DAAC, 2011).The root distribution in ISAM-Static is calculated based on the root depths at which plants have 50 % of their total root biomass and a dimensionless shape parameter for describing root profile (Schenk and Jackson, 2002).Since the static root distribution case assumes no temporal variation in root fraction in each soil layer, we use the average value of three observed corn root profiles (see Sect. 3) to calibrate the static root distribution case.The fixed canopy heights in the ISAM-StaticH experiment are assumed to be the maximum canopy height of specific vegetation type (H max ) from AmeriFlux data sets (Table A1).
In order to evaluate the performance of integrated effects of the dynamic crop growth processes implemented in this study (ISAM-Dynamic case) and the individual dynamic crop growth processes, we compare the Willmott indexes (dr d ) for carbon and energy fluxes based on individual five experiments discussed above with the estimated dr d for the ISAM-Dynamic case (Table 4).

Static vs. dynamic crop growth processes
The Willmott index values (dr d ) for daily mean GPP, Rn, SH and LH fluxes in ISAM-Dynamic case are higher than that in ISAM-Static case, and several are much closer to 1, except for no apparent improvement in dr d values for corn GPP and Rn fluxes at the Bondville site (Table 4).These results suggest that the implementation of dynamic crop growth scheme in ISAM significantly strengthens the ability of model to capture seasonal variability in measured carbon and energy fluxes for crops.No differences in dr d values for corn GPP and Rn fluxes at the Bondville site for ISAM-Dynamic and ISAM-Static experiments are due to that fact that processes considered in both experiments are unable to capture a crop lodging effect, as discussed in Sect.4.1.

Static vs. dynamic carbon allocation
Figures 1b, e, h, and k show that the estimated aboveground biomass for corn and soybean is in much better agreement with measurements for the ISAM-Dynamic case than for the ISAM-StaticC case.In addition, the ISAM-Dynamic case better captures the seasonal variability in leaf carbon mass, as indicated by LAI (Fig. 1a, d, g, j  biomass (Fig. 1h, k) than the ISAM-StaticC case.The improvements in estimated seasonal aboveground biomass, leaf and root carbon biomass for ISAM-Dynamic case are more for soybean than for corn at both sites.These results indicate that the dynamic carbon allocation scheme in the ISAM-Dynamic case is able to capture the response of carbon allocation to water, temperature, and light stresses, leading to a better simulation of aboveground total biomass and leaf carbon amount.With better simulated seasonal variability in carbon allocations, the dr d values for GPP, SH, and LH calculated based on ISAM-Dynamic case are generally closer to 1 than based on ISAM-StaticC case (Table 4), except for corn GPP at the Bondville site.No improvement in corn GPP at Bondville for ISAM-Dynamic is because the model is unable to capture the sharp reduction in GPP due to crop lodging with gusty wind, as discussed in Sect.4.1, even after accounting for the dynamic processes.Nevertheless, our results suggest that implementation of the dynamic carbon allocation parameterizations improves the model-estimated results for GPP, SH and LH fluxes, especially for soybean.

Static vs. dynamic LAI
Figure 1a, d, g, and j show that prescribed LAI usually underestimates LAI over the growing seasons at both the Mead and Bondville sites.In addition, prescribed LAI is not able to partition ground vegetation LAI and crop LAI, leading to wrong estimates of growing season length for the crop.The underestimation of the LAI over the growing season results in underestimation of the amount of solar radiation absorbed by the canopy, leading to underestimation of GPP and LH, but overestimation of SH.In contrast, the ISAM-Dynamic version of the model, which accounts for the dynamic green and brown LAI parameterizations, is able to capture observed seasonal variability in LAI (Fig. 1a, d, g, j).As a result of this, ISAM-Dynamic-based GPP, Rn, SH and LH fluxes for corn and soybean at both sites are in much better agreement with the observations than in the case of ISAM-StaticLAI, except for corn GPP and Rn at the Bondville site.The dr d values for ISAM-Dynamic are higher by 2-13 % for Rn, 3-41 % for GPP, 18-39 % for SH and 19-35 % for LH at both sites than for the ISAM-StaticLAI case (Table 4).The improvement for soybeans is usually larger than for corn.The smaller improvement for corn GPP and Rn in Bondville can be attributed to the fact that the ISAM-Dynamic and ISAM-Static cases are both unable to capture the effect of gusty wind on LAI.

Static vs. dynamic root distribution
In order to illustrate the importance of dynamic root characteristics, here we first compare the model estimated water uptake for the ISAM-Dynamic and ISAM-StaticR cases for 1980 at the Mead site.As discussed in Sect.3.1, we use 1980 NLDAS-2 climate forcing data (Mitchell et al., 2004) to drive the model.All other information, such as manage-ment seedling rate, planting time, etc., is taken from Newell and Wilhelm (1987).The ISAM-Dynamic case captures well the measured trends of root growth with soil depth during the growing season for corn at the Mead site; whereas the ISAM-StaticR case overestimates root density in shallow soil layers but underestimates in deep soil layers during the growing season (Fig. 5a-c).These differences in root characteristics for the two parameterization cases in turn influence the estimates of simulated soil water stress and root water uptake and hence transpiration (Fig. 5d-f).This is because transpiration is more sensitive to the moisture content of the densely rooted shallow soil layers than that in the remainder of the root zone (Feddes et al., 2001).On Julian day 174 when soil moisture is optimal during the initial vegetative stage, the calculated total amount of root water uptake for both cases is approximately the same (2.33 mm day −1 ), but there are substantial differences in the magnitude of the water uptake at different soil depths, specifically in the shallow soil layers.For the ISAM-StaticR case the maximum amount of water is extracted from shallow soil layers above 0.03 m, whereas for the ISAM-Dynamic case roots take water from the more moist deeper layers above 0.12 m (Fig. 5d).However, ISAM-StaticR parameterization overestimates the root density and water uptake in shallow layers (Fig. 5b-c), reducing the soil water available in the shallow soil as the growing season progresses.This results in an earlier and more intense start of soil moisture stress and lower actual transpiration in the ISAM-StaticR case than that in the ISAM-Dynamic case (Fig. 5e-f).In order to illustrate the importance of the dynamic root distribution scheme for seasonal variability in crop transpiration, here we compare the simulated transpiration results for the 2001 corn-growing season at Mead.In Fig. 6, the transpiration is higher for ISAM-Dynamic than for the ISAM-StaticR case during the growing season, and the transpiration differences between the two cases gradually increase, especially during the summer, when low summer precipitation cannot effectively compensate soil water depletion in shallow layers (not shown here).The increased water uptake from deeper and moister root zones in the ISAM-Dynamic case mitigates the intensity of water stress during the growing season by about 60 % of that of the ISAM-StaticR case and improves the simulations of soil water uptake when soil water in the upper soil layers is exhausted during the growing season.
In order to evaluate the validity of the dynamic root parameterization scheme, we compare model results for total transpiration, latent heat flux, and GPP during the 2001-2004 growing seasons with corn and soybean rotations at the Mead and Bondville sites.The ISAM results suggest that ISAM-Dynamic parameterization-estimated plant water transpirations during the 2001-2003 growing season are about 28-34 % higher than ISAM-StaticR (Fig. 7a).However, there is no apparent difference in plant water transpiration between the ISAM-StaticR and ISAM-Dynamic cases over the 2004 growing season at both sites (Fig. 7a).Both sites experience moist weather conditions during the summer of 2004.For example, the accumulated precipitation rates at the Mead and Bondville sites for the period of June to July 2004 are about 91 and 63 % higher than the average for the same time period over 2001-2003.Therefore, soybean experiences no water stress conditions at either site in 2004, and the estimated water transpiration fluxes for ISAM-StaticR and ISAM-Dynamic cases are approximately the same (Fig. 7a).This also results in similar values for estimated GPP and LH fluxes for both cases (Fig. 7b-c).The increased transpiration in the ISAM-Dynamic case, relative to the ISAM-StaticR case, mitigates the water stress effect on catalytic capacity of RuBisCO (V c max25 ) and stomatal conductance.This results in a 13-61 % increase in GPP and 12-27 % increase in LH at the Mead site, and a 26-41 % increase in GPP and 13-21 % increase in LH at the Bondville site for the ISAM-Dynamic case relative to the ISAM-StaticR case (Fig. 7b-c).The increased values for GPP and LH for the ISAM-Dynamic case are in much better agreement with observations (Fig. 7b-c) than for the ISAM-StaticR case.Moreover, the dr d and dr h values for GPP and LH (Table 4) are much closer to 1 in the ISAM-Dynamic case than those in the ISAM-StaticR case for most of the cases.One particular exceptional case is corn GPP at the Bondville site.As discussed in Sect.4.1, ISAM-Dynamic is not able to capture a sharp reduction in corn GPP during this year, and thus overestimates corn GPP at the Bondville site.The ISAM-Dynamic case mitigates downscaled effects of water stress on GPP and thus further overestimates GPP (Fig. 7b-c).The increase in dr d values for GPP and LH for the ISAM-Dynamic case as compared to the ISAM-StaticR case (Table 4) suggests that ISAM-Dynamic much better captures the seasonal pattern of carbon, energy and water fluxes than the ISAM-StaticR.Specifically, it much better captures the apparent increase in the values of dr d for GPP and LH at the Mead site (Table 4), where crops endured water stress conditions during the 2001-2003 growing season, indicating the importance of dynamic carbon allocation and root distribution mechanisms in the calculations of carbon, energy, and water fluxes under water stress conditions.

Static vs. dynamic canopy height
Table 4 shows that dr d values have small differences between ISAM-StaticH and ISAM-Dynamic cases, relative to comparisons discussed above, indicating that the implementation of dynamic canopy height simulation does not apparently improve the carbon and energy fluxes for these crops.This is perhaps due to the fact that there is no large seasonal variability in canopy height for corn and soybean.Thus, replacing prescribed canopy height with seasonally variable canopy height does not significantly change the atmospheric turbulence above the crop canopy or the carbon and energy fluxes.

Conclusions
We have implemented dynamic crop growth processes into a land surface model, ISAM.These dynamic crop growth processes include specific phenology development for corn and soybeans, dynamic carbon allocation, vegetation structure and root distribution, as well as different removal rates for fresh and old standing brown leaves and the effects of leaf age on CO 2 compensation points for C3 crops.The C3 and C4 crop growth processes in the model are calibrated with half-hourly/hourly data for LAI, biomass, and carbon, water, and energy fluxes measured for corn-soybean rotation systems at the Mead, Nebraska AmeriFlux site, and the model is evaluated for the same variables using the data from another AmeriFlux site in Bondville, Illinois.The calibrated and evaluated ISAM is able to capture the diurnal and growing season patterns of carbon assimilation, and water and energy fluxes for corn (C4 crop) and soybean (C3 crop) at these two sites.Specifically, the calculated GPP, Rn, and LH fluxes compared well with observations, but the model is unable to capture the variation in SH flux for corn and soybean at both sites as discussed in Sect.4.1.
The model's dynamic carbon allocation parameterization, dynamic LAI and phenology development, dynamic canopy height, and dynamic root distribution capture the measured seasonal patterns of vegetation structures well, in particular changes in LAI and the vertical distribution of roots in soil.The improvement in vegetation structure simulation better captures the seasonal variability in carbon and energy fluxes at both sites, relative to the static simulations of vegetation structure.With dynamic carbon allocation and root distribution schemes, the improved crop water transpiration and soil water stress significantly improve modeled GPP and LH, especially during dry periods.The percent differences between the estimated fluxes based on the ISAM-Dynamic and ISAM-Static cases for LH are 12-27 % and for GPP 13-62 % at the Mead and Bondville sites.These results indicate the importance of considering dynamic allocation and root distribution processes in land surface models to simulate the carbon, water, and energy fluxes accurately, especially during dry periods.The incorporation of the effect of leaf age on Overall, the implementation of dynamic crop growth processes into the ISAM allows us to study the feedbacks between carbon, water, and energy fluxes and environmental variables.The extended ISAM can be applied on a large scale, where carbon, water, and energy fluxes can vary with spatial variations in environmental variables.Furthermore, the implementation of dynamic crop growth processes is particularly important for understanding the interactions between the land surface processes and climate change.For example, with a dynamic carbon allocation and vegetation structure scheme, the extended ISAM can simulate the adaptation of vegetation to climate change.These adaptations include adjustments in behavior, morphology, and physiology.(Hendry et al., 2008).Due to the lack of this scheme, current land surface models possibly overestimate climate effects on terrestrial ecosystems.In addition, with dynamic root distribution schemes, the model is able to simulate the root response to soil water availability in each soil layer and accounts for redistribution of soil water by root systems.The implementation of this scheme can improve the simulation of vegetation transpiration under dry conditions by extracting water from deeper and moister soil layers.Jasechko et al. ( 2013) have found that vegetation transpiration represents 80-90 % of terrestrial water fluxes but is underestimated by the current land surface models.With dynamic treatment of carbon allocation and root distribution, it is expected that the extended ISAM can correct currently underestimated global vegetation transpiration and overestimated soil and canopy evaporation.
In future studies, the model will be applied to assess the interaction between crop growth and climate change.Since we have developed a flexible process-based crop growth modeling framework, it can also be applied to simulate not only other food crops, such as wheat, but also energy crops, such as Miscanthus and switchgrass.

Carbon
For carbon allocation during initial and post-reproductive period Eq.(A26) allocation when NPP > 0  Root depth and Root fraction in each soil layer (j ) Eq. (A43) distribution r j = f j j = 1 f j − f j −1 2 ≤ j ≤ L max

Fig. 1 .
Fig. 1.Measured and model simulated LAI, aboveground biomass, root biomass, canopy height under corn and soybean rotation at the Mead and Bondville AmeriFlux sites.Measured data for corn root biomass are available for the 2001 growing season and for soybean over the 2002 and 2004 growing seasons at the Bondville, IL, site.The top and the bottom panels for corn column are for 2001 and 2003 growing seasons, whereas for soybean column are for 2002 and 2004 growing seasons.

44
Figure 2. Measured and model simulated mean diurnal variations in gross primary productivity (GPP), net radiation (Rn) at the canopy top, sensible heat (SH), and latent heat (LH) fluxes.The diurnal cycle of each flux shown here for corn and soybean is the mean diurnal cycle over two growing seasons.For corn, the diurnal cycle is averaged over 2001 and 2003 growing seasons, whereas for soybean it is averaged over 2002 and 2004 growing seasons.The error bars indicate ± 1 standard deviation (SD) of variation for hourly/half hourly values over the two growing seasons.

Fig. 2 .
Fig. 2. Measured and model simulated mean diurnal variations in gross primary productivity (GPP), net radiation (Rn) at the canopy top, sensible heat (SH), and latent heat (LH) fluxes.The diurnal cycle of each flux shown here for corn and soybean is the mean diurnal cycle over two growing seasons.For corn, the diurnal cycle is averaged over 2001 and 2003 growing seasons, whereas for soybean it is averaged over 2002 and 2004 growing seasons.The error bars indicate ±1 standard deviation (SD) of variation for hourly/half-hourly values over the two growing seasons.

Figure 3 .
Figure 3. Measured and model simulated daily mean gross primary productivity (GPP), net radiation (Rn) top of the canopy, sensible heat (SH), latent heat (LH), and soil water (SW) under corn and soybean rotation at Mead and Bondville over 2001-2004 growing seasons.Flux values for individual sites are represented by a set of two figures.For corn, the top panel figure shows the flux values for the 2001 growing season and the bottom panel for the 2003, whereas for soybean the top panel figure shows the flux values for the 2002 growing season and the bottom panel for the 2004.

Fig. 3 . 1160 Fig. 4 .
Fig. 3. Measured and model simulated daily mean gross primary productivity (GPP), net radiation (Rn) top of the canopy, sensible heat (SH), latent heat (LH), and soil water (SW) under corn and soybean rotation at Mead and Bondville over 2001-2004 growing seasons.Flux values for individual sites are represented by a set of two figures.For corn, the top panel figure shows the flux values for the 2001 growing season and the bottom panel for the 2003, whereas for soybean the top panel figure shows the flux values for the 2002 growing season and the bottom panel for the 2004.
Figures 1b, e, h, and kshow that the estimated aboveground biomass for corn and soybean is in much better agreement with measurements for the ISAM-Dynamic case than for the ISAM-StaticC case.In addition, the ISAM-Dynamic case better captures the seasonal variability in leaf carbon mass, as indicated byLAI (Fig. 1a, d, g, j), and the root carbon

Figure 5 .Fig. 5 .
Figure 5.Comparison of modeled and measured corn root density (a-c) and water uptake (e-f) 1161 profiles (0-2 m) for three different days during the growing season at the Mead site.The data 1162 and model results for the dynamics (ISAM-Dynamic) and Static (ISAM-StatcR) cases are plotted 1163 for 1980.1164 1165

1166Figure 6 .Fig. 6 .
Fig. 6.Model estimated daily water uptake for the ISAM-Dynamic and ISAM-StaticR cases during the 2001 corn growing season at the Mead site.

Table 4 .
The Willmott index (dr d ) to quantify the degree to which observed daily mean GPP and energy fluxes are captured by the model for corn and soybean at the Mead and Bondville sites.The n is the number of observation at the daily step.

Fig. 7 .
Fig. 7. Measured and simulated total GPP, transpiration and latent heat fluxes (LH) from 2001-2004 growing season under corn-soybean rotation at the Mead and Bondville sites.The odd years 2001 and 2003 are the corn planting years, whereas the even years 2002 and 2004 are the soybean planting years.

Table 1 .
Calibrated processes and parameters and their original and updated values.The two data values in original and calibrated columns are for corn and soybean, respectively.

Table 2 .
The Willmott index to quantify the degree to which observed GPP, energy and water fluxes are captured by the model for corn and soybean at the Mead and Bondville site.The dr h is Willmott index for hourly observed data and model results, and dr d is index for daily mean observed data and model results.The N h is the number of observation samples at the hourly time step, and the N d is the number of observation at the daily time step.
(Winstanley, 2003)d values for corn at the Bondville site are lower than the dr h values, indicating that the model estimates are less consistent with the measured values during a certain period of the growing season.Figure3ksuggests that the model fails to capture a sharp reduction in GPP during the initial reproductive period of 2003 (between Julian day 202 and 215).The reason for the sharp reduction in observed GPP is unknown, but the Illinois water and climate summary on July 2003 reports widespread crop lodging due to gusty wind during this period in central Illinois(Winstanley, 2003).The weather report at the nearest weather station (40.03 • N, 88.28 • W)(Climate Champaign,  2003)also suggests that the area received a thunderstorm with wind gusts over 30 miles h −1 on Julian day 202 of 2003 and wind gusts (> 13.5 m s −1 ) between Julian day 203 and 208 of 2003.The high wind gusts might have induced crop lodging and hence reduced the GPP.The model is unable to capture this information, because the model does not currently account for the effect of extreme wind gusts on crops.

Table 3 .
Description of the model experiments performed to evaluate the effects of different dynamic processes on model results.

Y. Song et al.: Implementation of dynamic crop growth processes into a land surface model Appendix ATable A1 .
Variables and parameters that appear in the model equations.

Table A2 .
ISAM equations in this study.

Table A2 .
Continued.= GPP i − (R m_leaf i + R m_stem i + R m_root i + R m_grain i ) − R g i R m_leaf i = k leaf * C g_leaf i +C d_leaf i CN leaf * g(T leaf ) R m_stem i = k stem * C stem CN stem * g(T leaf ) R m_root i = k root * C root CN root * g (T soil ) R m_grain i = k grain * C grain CN grain * g(T leaf ) 10_above = 3.22 − 0.046 * (T leaf − 273.16)Q 10_below = 3.22 − 0.046 * (T soil − 273.16) = max 0, 0.25 * GPP i − R m leaf i − R m stem i − R m root i − R m grain i R g leaf i = R g i * C g leaf i +C d leaf i C g leaf i +C d_leaf i +C stem i +C root i +C grain i R g stem i = R g i * C stem i C g leaf i +C d_leaf i +C stem i +C root i +C grain i R g root i = R g i * C root i C g leaf i +C d_leaf i +C stem i +C root i +C grain i R g grain i = R g i * C grain i C g leaf i +C d_leaf i +C stem i +C root i +C grain i