Development of the DayCent-Photo model and integration of variable photosynthetic capacity

Abstract

We integrated a photosynthetic sub-model into the daily Century model (DayCent) to improve the estimations of seasonal changes in carbon fluxes at the Niwot Ridge LTER site and the Harvard forest LTER site (DayCent-Photo). The photosynthetic sub-model, adapted from the SIPNET/PnET family of models, includes solar radiation and vapor pressure deficit controls on production, as well as temperature and water stress terms. A key feature we added to the base photosynthetic equations is the addition of a variable maximum net photosynthetic rate (Amax). We optimized the parameters controlling photosynthesis using a variation of the Metropolis-Hastings algorithm along with data-assimilation techniques. The model was optimized and validated against observed net ecosystem exchange (NEE) and estimated gross primary production (GPP) and ecosystem respiration (RESP) values for AmeriFlux sites at Niwot Ridge and Harvard forest. The inclusion of a variable Amax rate greatly improved model performance (NEE RMSE= 0.63 gC·m⁻², AIC= 2099) versus a version with a single Amax parameter (NEE RMSE= 0.74 gC·m⁻², AIC= 3724). DayCent-Photo was able to capture the inter-annual and seasonal flux patterns for NEE, GPP, ecosystem respiration (RESP), and daily actual evapotranspiration (AET), but tended to overestimate yearly NEE uptake. The DayCent-Photo model has been successfully set up to simulate daily NEE, GPP, RESP, and AET for deciduous forest, conifer forests, and grassland systems in the US using AmeriFlux data sets and has recently been improved to include the impact of UV radiation surface litter decay (DayCent-UV). The simulated influence of a variable Amax rate suggests a need for further studies on the process controls affecting the seasonal photosynthetic rates. The results for all of the forest and grassland sites show that maximum Amax values occurs early during the growing period and taper off toward the end of the growing season.

Keywords

DayCent-Photo model / seasonal maximum net photosynthetic rate / net ecosystem exchange / gross primary production / UV radiation

Introduction

Forest ecosystems exchange vast amounts of CO₂ to and from the atmosphere (^{Schimel, 1995}), serving as a net sink (^{Luyssaert et al., 2008}, ²⁰¹⁰), and playing an important role in curbing the increase of atmospheric CO₂ concentrations. In forest ecosystems, the net exchange of carbon between atmosphere and live biomass is strongly regulated by the seasonality of photosynthesis (^{Huxman et al., 2003}; ^{Monson et al., 2005}). However, the seasonality is poorly represented in terrestrial biosphere models (^{Sacks et al., 2006; Richardson et al., 2013}). The AmeriFlux net carbon exchange studies provide a data set to measure and understand the factors affecting net ecosystem exchange (NEE) in different biomes across the globe. This study utilizes the long term forest ecosystem observations of the Niwot Ridge AmeriFlux site to explore the role of variable photosynthetic rates on seasonal carbon exchange for conifer systems and shows how variable maximum photosynthesis rates impact ecosystem dynamics for deciduous forest and grassland systems using the photosynthetic version of the DayCent ecosystem model (DayCent-Photo) (^{Savage et al., 2013}).

The seasonality of photosynthesis results from changes in photosynthetic capacity, and affects gross primary production (GPP) and ultimately net ecosystem exchange. Aspects including the variation of stem xylem hydraulic conductivity, enzyme activity, mesophyll conductance, stomatal conductance, electron transport efficiency, quantum efficiency, and general leaf design principles have been found to relate to the seasonal dynamics of photosynthesis (^{Granda et al., 2014}; ^{Urban et al., 2017}; ^{Zhang et al., 2017}; ^{Martinez and Fridley, 2018}). The seasonal changes in photosynthetic capacity are well-documented in forest systems (^{Bourdeau, 1959}; ^{Helms, 1965}; ^{Mohren and van de Veen, 1995}; ^{Marshall et al., 2001}; ^{Huxman et al., 2003}; ^{Wang and Barrett, 2003}; ^{McGarvey et al., 2004}; ^{Weiskittel et al., 2006}). In an evergreen forest, photosynthetic capacity increases rapidly in the spring warm up, allowing the forest to assimilate a large portion of GPP for the year during the early snow melt period (^{Monson et al., 2005}). Photosynthetic capacity then generally declines in the dry or cool season in most (woody) species (^{Zhang et al., 2014}, ²⁰¹⁷). The seasonality of photosynthesis varies among species. ^{Guan et al. (2014)} found that the net photosynthetic rate peaked in summer for two species (Pinus massoniana and Castanopsis sclerophylla) but in autumn for two other species (Lithocarpus glaber and Cyclobalanopsis glauca). ^{Urban et al. (2017)} found that photosynthetic capacity peaked in August and September for Norway spruce (i.e., Picea abies). ^{Li et al. (2018)} reported that for Ammopiptanthus mongolicus, photosynthetic capacity increased from January, reached the highest value in July, and greatly declined in September. ^{Martinez and Fridley (2018)} found that photosynthetic capacity peaked in spring and declined through summer and autumn for common species in central New York (i.e., Acer negundo, Liriodendron tulipifera, Acer saccharum, Tsuga Canadensis, and Pinus strobus). The seasonality of photosynthesis and other ecological processes is usually examined with short-term experiments or observations at small plots, and projecting to larger temporal and spatial scales involves mathematical models such as the DayCent model (^{Rastetter et al., 2003}). Incorporating the seasonality of photosynthetic capacity to terrestrial biosphere models will allow the models to better simulate the inter-annual dynamics in NEE and improve carbon flux estimates at the seasonal and annual level (^{Wang et al., 2007}; ^{Tang et al., 2013}).

The original version of the DayCent model (^{Parton et al., 2001}) was developed to model soil trace gas fluxes such as N₂O, NO_x, CH₄, and CO₂, and plant production for forest, cropland, grassland, and savanna ecosystems. A version specific to forests was developed to better represent the ecosystem dynamics for forest systems (ForCent) (^{Parton et al., 2010}). All DayCent versions include complex controls of soil nutrients, soil water and temperature, and atmospheric CO₂ levels on plant growth (^{Parton et al., 2010}). The DayCent model have been used to simulate long-term ecological processes at forest sites such as the Harvard Forest site in Massachusetts and the Niwot Ridge site in Colorado. ^{Frey et al. (2013)} used the DayCent model to simulate the effects of soil warming on soil C storage at the Harvard forest site in response to altering the microbial efficiency parameter in the model. ^{Savage et al. (2013)} used the DayCent model to simulate the long-term impacts of warming and N deposition on the distribution of forest C stocks at the Harvard forest site. ^{Hartman et al. (2014)} used the DayCent model to investigate the combined effects of climate, atmospheric N deposition, and CO₂ on many ecosystem variables at seven forested montane including the Niwot Ridge site.

The DayCent model was updated to include a photosynthetic sub-model from the Simple Photosynthesis and Evapotranspiration (SIPNET) model (^{Braswell et al., 2005}; ^{Sacks et al., 2006}) to better simulate ecosystem dynamics (^{Savage et al., 2013}), and modified to incorporate the seasonally variable photosynthetic capacity in the calculation of GPP (DayCent-Photo). This paper will present a detailed description of the DayCent-Photo model and its validation against eddy covariance data collected at Niwot Ridge AmeriFlux site. DayCent-Photo has recently been revised to include the impact of ultraviolet-B (UV-B) radiation on grassland and forest ecosystems (DayCent-UV) and has been successfully used to simulate ecosystem dynamics for the short grass steppe ecosystem in Eastern Colorado (^{Chen et al., 2016}, ²⁰¹⁷). The data collected at the Niwot Ridge LTER high-altitude forest site include extended series of carbon and water flux measurements, and have been successfully used with model-data fusion techniques to provide well constrained estimates of system level fluxes (^{Braswell et al., 2005}; ^{Sacks et al., 2006}, ²⁰⁰⁷; ^{Moore et al., 2008}; ^{Zobitz et al., 2008}). We parameterized DayCent-Photo with a model-data synthesis technique that is a variation of the Markov Chain Monte Carlo (MCMC) algorithm, and evaluated model performance with a combined likelihood given daily fluxes observed at the Niwot Ridge site. We will compare two version of the model, one with variable maximum photosynthetic capacity (Amax), and one with a fixed Amax. Since we don’t have a variable Amax function, optimizing the parameters controlling photosynthesis, using model-data fusion techniques, will create a synthetic stand-in for that function. We will then do a model comparison to see if variable Amax could be a significant factor in control production.

Methods

Here we describe the DayCent-Photo model which includes the SIPNET photosynthesis model, the NIWOT ridge high elevation conifer site, the eddy covariance data used to calibrate and test the model performance, and the procedure used to determine the optimal parameters for the SIPNET sub-model. We used the observed daily value of NEE, GPP, and actual evapotranspiration (AET) from the eddy covariance site at NIWOT ridge to calibrate and test the model.

Niwot Ridge site and data description

Data used were from the Niwot Ridge AmeriFlux site (US-NR2) (40.0329°N Latitude; 105.5464°W Longitude), a high-elevation (3030 m), subalpine forest on the western side of the Continental Divide. The site is a part of the Roosevelt National Forest, approximately 50 miles west of Boulder, Colorado in the Rocky Mountain Range, with mean annual temperature of 4°C and mean annual precipitation of 800 mm of which ~65% falling as snow (^{Monson et al., 2002}). The dominant tree species in this area are subalpine fir (Abies lasiocarpa), Engelmann spruce (Picea engelmannii), and lodgepole pine (Pinus contorta). The forest stand at the site is approximately 100 years old and is a natural secondary growth, established following clear-cut logging from 1900 to 1910. The site’s understory is sparse, with 25% average cover, and is composed of seedlings, Whortleberry (Vaccinium myrtillus L.), moss, and lichens. The soil at the site is sandy due to rocky granite composition, with a surface layer of ~6 cm of organic material. Estimates for above ground biomass of carbon range from 11,500 to 24,000 gC·m^‒2 (^{Blanken et al., 2010}), and the site is a carbon sink with a cumulative annual rate ranging from 60‒80 gC·m^‒2 (^{Monson et al., 2002}; ^{Tang et al., 2012}).

The data at the site were collected by eddy covariance method (^{Baldocchi et al., 1988; Baldocchi, 2003}) and retrieved from the AmeriFlux website (^{Blanken et al., 2010}). Briefly, the eddy covariance method measured CO₂ concentration above the canopy with a sonic anemometer and an infrared gas analyzer, converted the raw data to half hourly fluxes of NEE applying various correction measures. The data were gap-filled using the Marginal Distribution Sampling method (^{Papale and Valentini, 2003}; ^{Reichstein et al., 2005}) and u-star filtered, to account for periods of low turbulence in the air mass that affect flux calibrations. Estimates for ecosystem respiration (RESP) were calculated by the short-term temperature response of nighttime NEE fluxes (^{Reichstein et al., 2005}), and GPP estimates were calculated as GPP=RESP‒NEE. The data were in half hour flux rates, and were aggregated into daily values. A detailed description of data collection is available elsewhere (^{Monson et al., 2002}; ^{Turnipseed et al., 2002}, ²⁰⁰³, ²⁰⁰⁴). The climate variables were air temperature, precipitation, global radiation, and vapor pressure deficit while the flux variables included GPP, NEE, RESP, and AET.

Model structure

The DayCent-Photo model is an ecosystem model (Fig. 1) capable of simulating plant growth, soil carbon processes, daily NEE, GPP, and AET fluxes, soil nutrient cycling (N, P, S, and K) and trace gas fluxes (N₂O, NO_x, N₂, and CH₄), and is set up to simulate the impact of management practices on natural and managed grasslands, forests, cropland, and savannas ecosystems (^{Parton et al., 1998, 2001}; ^{Del Grosso et al., 2000}). Plant growth is simulated through the processes of assimilation, allocation of carbon to different plant parts, and exchange of water, carbon and other nutrients. The processes occur in various compartments of the system, with each compartment representing a different plant and soil functions (roots, wood, leaves, litter layer, soil, etc.) The earlier version of the ForCent model (^{Parton et al., 2010}) estimates GPP as a scalar value twice that of the calculated NPP, and NPP in turn is a function of genetic potential, soil temperature, available water, solar radiation, and available nutrients (^{Parton et al., 2010}).

The DayCent-Photo model uses the SIPNET with simplified Farquar photosynthesis to calculate daily GPP as a function of maximum photosynthetic assimilation rate (GPP_max) that is the unstressed (not including climatic factors) net instantaneous CO₂ assimilation rate per unit leaf mass:

(1){{\displaystyle \mathrm{G}\mathrm{P}\mathrm{P}}}_{\max }=A\mathit{max}*(A\mathit{max}Frac+BaseFol\mathit{Re}spFrac),

where Amax is the maximum net CO₂ assimilation rate per unit leaf mass, AmaxFrac is the average daily maximum photosynthesis as a fraction of Amax, and BaseFolRespFrac is the foliar respiration as a fraction of Amax.

The GPP_max value is then scaled by the four terms that represent climatic conditions with values between zero and one: temperature (D_temp), vapor pressure deficit (D_vpd), solar radiation (D_light), and a water stress term (D_water) such that

(2)\mathrm{G}\mathrm{P}\mathrm{P}={{\displaystyle \mathrm{G}\mathrm{P}\mathrm{P}}}_{\max }*{{\displaystyle D}}_{\mathrm{t}\mathrm{e}\mathrm{m}\mathrm{p}}*{{\displaystyle D}}_{\mathrm{l}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}*{{\displaystyle D}}_{\mathrm{w}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{r}}.

This equation is from the SIPNET model (^{Monson et al., 2002}; ^{Braswell et al., 2005}; ^{Sacks et al., 2006}, ²⁰⁰⁷) and the SIPNET model is based on the PnET model (^{Aber and Federer, 1992}) designed to be used with eddy covariance data in a model-data fusion approach for estimating parameter values. Further descriptions and equations for these scaling terms can be found in ^{Braswell et al. (2005)}. We modified the calculation for D_vpd to improve the effect of stomatal conductance due to lack of air moisture (VPD). The values for D_vpd is calculated as an exponential decline in efficiency (^{Drake et al., 2010}):

(3){{\displaystyle D}}_{\mathrm{v}\mathrm{p}\mathrm{d}}=DVPDSlope*e^{DVPDExp*VPD},

where DVPDExp and DVPDSlope are the shape parameters controlling the VPD photosynthesis relationship.

The original version of the SIPNET model included a single Amax throughout all seasons. We changed the SIPNET model so that seasonal changes in Amax are represented using a linear Amax response function as observed by ^{Marshall et al. (2001)}in Douglas-fir, Ponderosa pine, and Western white pine of northern Idaho. The response function scales the value of Amax from the start of one growth period to the end of the growing period (Fig. 2) as defined by the GrowthDays parameters (the number of days since the start of the growing season). The function allows the Amax value to be scaled three times during the growing season. Figure 2 shows the seasonal changes in the Amax scalar derived from the model optimization process. The figure also shows the seasonal changes in the Amax scalar for the Harvard forest site, and the SGS short grass steppe site (^{Savage et al., 2013}; ^{Chen et al., 2017}). The optimal values for the Harvard forest and SGS sites were determined using eddy covariance data sets.

Statistical analysis and parameter optimization

We used an optimization technique to estimate the parameter values for both the SIPNET photosynthesis routine and for the variable Amax function, optimized against observations of GPP and NEE. The additional use of GPP and NEE for optimization greatly improves the model results for all flux estimates and keeps the parameters controlling photosynthesis from over-fitting the data and exceeding the literature-defined parameter boundaries. We chose not to use AET as an optimization parameter because changes in photosynthetic parameters do not affect the value of AET in our version of the photosynthetic sub-model. The parameter optimization scheme used is a variation of the Metropolis-Hastings algorithm, with a simulated annealing algorithm (^{Hurtt and Armstrong, 1996}) to help prevent parameter values from settling into local optima. The algorithm performs a quasi-random walk through the parameter space, via the Metropolis-Hastings algorithm, to find a parameter set that minimizes model error and maximizes the likelihood function (L) (^{Braswell et al., 2005}; ^{Sacks et al., 2006, 2007}; ^{Moore et al., 2008}). The algorithm produces a maximum likelihood estimation (MLE), which is the parameter vector that most strongly supports the data (^{Fisher, 1932}).

This model optimization method based on maximum likelihood estimation allows for several candidate models (or model parameterizations in this case) to be evaluated and ranked based on the support by the data (^{Johnson and Omland, 2004}). The likelihood function compares model hypothesis by calculating confidence bounds of parameters, which we cannot do with sum of squares (^{Hilborn and Mangel, 1997}).

We sum the calculated probabilities for multiple parameters by using log-likelihood for a total likelihood or Ltotal. By convention we are using a negative log-likelihood, as seen in the following equation, since we are maximizing the likelihood function. Future references to log-likelihood should be considered a negative log-likelihood. The best parameter set produces a result that maximizes log-likelihood which is computed as

(4)l=[n\hspace{0.17em}\log \hspace{0.17em}\sigma ]-\left[{\displaystyle \frac{\mathrm{1}}{{\mathrm{2}\sigma }^{\mathrm{2}}}}{{\displaystyle \Sigma }}_{i=\mathrm{1}}^n{({{\displaystyle x}}_i-{{\displaystyle \mu }}_i)}^{\mathrm{2}}\right],

where n is the number of observations, s is the standard deviation for the data, {{\displaystyle x}}_i is an observed data point, and {{\displaystyle \mu }}_i is the corresponding modeled point. This form of the log-likelihood equation is easy to work with because the number of observations n and the sum of squares of the residuals is readily available, and the standard deviation can be estimated from the sum of squares. Constants that do not affect the likelihood comparison were removed from the equation. Standard deviation, s, was estimated because the actual error for NEE and GPP data is unknown. We estimated the standard deviation, s_e, according to the method in ^{Braswell et al. (2005)}:

(5){{\displaystyle \sigma }}_e=\sqrt{{\displaystyle \frac{\mathrm{1}}{n}}{{\displaystyle \Sigma }}_{i=\mathrm{1}}^n{({{\displaystyle x}}_i-{{\displaystyle \mu }}_i)}^{\mathrm{2}}}.

The optimization process proceeds in iterative sequences. For each run (a 2000 years’ simulation), the optimization process randomly selects a parameter, changes it by a random amount, and compares the resulting model output against observed NEE and estimated GPP data. If the new candidate parameter set increases log-likelihood, that set is then accepted. If the new candidate set does not have an improved likelihood, it may still be accepted with a probability equal to the difference of the log-likelihoods from the old to the new parameter sets to ensure that the global optimum is found. The likelihood calculation is an iterative Bayesian calculation sampled from a posterior distribution of the model error that is assumed to have a Gaussian probability distribution. After an adequate number of runs, the sampled space should represent the posterior probability distribution of parameters, and the Metropolis-Hastings algorithm guarantees convergence to a stationary point in the parameter space. Parameter ranges and initial starting values that are directly equivalent to the parameters of the SIPNET model are used from the previous work of ^{Moore et al. (2008)}, see Table 1. A near optimal parameter set, suitable for a model comparison, was retrieved with the methodology set forth in this paper. Extensive plant production, soil carbon and nitrogen, and biomass data from the Niwot Ridge LTER site (http://niwot.colorado.edu/) were used to calibrate the ability of the model to correctly simulate observed soil carbon and nitrogen, plant production and live plant biomass (^{Del Grosso et al., 2008}) levels for the Niwot Ridge site. R² and Root Mean Squared Error (RMSE) were used as the statistical method to validate the results. R² measures the proportion of variability in Y explained by the regression model. It is the square of the correlation r. RMSE measures the standard deviation of the residuals (the spread of the points about the fitted regression line).

Results

Optimization of photosynthesis parameters

We compared the ^{Moore et al. (2008)}best parameter estimates for the SIPNET photosynthesis model (Table 1) with the optimum SIPNET photosynthesis parameters in DayCent-Photo for the single and variable Amax runs (Table 2). As expected, these were significantly different from the original starting parameter used from ^{Moore et al. (2008)}. We expected to see differences in the parameterizations because the DayCent-Photo was optimized on two carbon flux values (NEE and GPP), whereas the SIPNET model was optimized on NEE and AET. Values of the best parameterizations fell well within the allowable range for each parameter for both DayCent-Photo variable Amax and single Amax versions and were generally quite similar for both versions of the model. Figure 2 shows the seasonal changes in the Amax scalar derived from the optimization of the photosynthesis parameters for the Niwot ridge site and three other sites from previous studies (^{Savage et al., 2013}; ^{Chen et al., 2016}, ²⁰¹⁷). The seasonal pattern of the Amax scalar was similar among DayCent-Photo (or DayCent-UV) adjusted for Niwot Ridge, Harvard forest (^{Savage et al., 2013}), short grass steppe of Colorado (CPER) (^{Chen et al., 2016}, ²⁰¹⁷), and tall grass prairie at Konza Biological Station in northern Kansas. The results show a rapid increase in the Amax scalar during the spring (April 1 to May 30), slow decrease during the summer, and a rapid decrease at the end of the growing season.

Simulation of observed seasonal patterns

Variable Amax improved the ability of DayCent-Photo to simulate observed seasonal patterns in NEE and GPP at Niwot Ridge (Table 3). Variable Amax was favored over single Amax by all statistical indices used (Table 3). DayCent-Photo with variable Amax had substantially lower RMSE for both NEE and GPP with the variable Amax version (0.63 vs. 0.74 gC·m^‒2 for NEE and 0.88 vs. 1.02 gC·m^‒2 for GPP) while the higher R² values for NEE and GPP. The AIC value for the variable Amax model (2099.0) was significantly better than for the single Amax model (3723.9) (Table 3).

With variable Amax, DayCent-Photo simulated the phenology of NEE and GPP reasonably well at the Niwot Ridge site, whereas the single Amax version greatly underestimated the spring uptake (Fig. 3). The comparison of observed and simulated mean weekly NEE showed a strong carbon uptake (downward pattern) during the early spring snow melt (weeks 16 to 20-May) for the observed NEE, followed by a reduction of the uptake in mid-summer as heterotrophic respiration increase, and another smaller increase in the latter summer period that coincides with summer monsoonal moisture common in this study area. Examination of NEE residuals between weeks 13‒25 showed over twice the cumulative error for the single Amax model (7.6 gC·m^‒2) compared to the variable Amax version (3.4 gC·m^‒2). Residual errors for the remainder of the growing season on were similar at 2.5 and 2.9 gC·m^‒2, respectively. The observations for GPP (Fig. 3(b)) showed a rapid increase during the spring (weeks 14 to 20), a gradual decline during the summer, and a rapid decline in the fall. The variable Amax version again simulated this phenology reasonably well, whereas the single Amax version underestimated the spring increase and did not show the summer decline. These results suggest that variable Amax improved the seasonal predictions of NEE by correctly simulating the rapid increase in GPP during the spring (weeks 14 to 20—April 15 to May 30). Another factor which contributed to the better results with the variable Amax version is the optimal temperature curve for the variable Amax version which shows that the variable Amax version temperature curve starts much earlier and at a higher scalar value (not shown). The net effect on GPP is a quicker and larger response to early season temperature increases for the variable Amax version.

DayCent-Photo model validation

Model validation was done with GPP and NEE data withheld from the optimization period from 2006 to 2008, except estimates of GPP were not available for 2008. Validation runs for the DayCent-Photo model were done using the best parameter set from the optimization process. A comparison for the optimization and validation period showed similar R² values for GPP and NEE estimates. The similar NEE and GPP results obtained (Fig. 4) for the optimization and validation period strongly support the addition of not only the photosynthetic sub-model, but also the inclusion of the seasonally variable Amax.

Figure 5 shows the comparison of observed vs. simulated daily NEE, GPP, and RESP for the Niwot Ridge site for 2005 using the variable Amax version of DayCent-Photo. The results for GPP and RESP (Fig. 5) show that the model simulated observed daily value quite well (model vs. observed data R²>0.88). The model simulated correct seasonal pattern in GPP with maximum values of GPP in the spring, rapid increases in GPP from April to May, and rapid decreases in GPP during the September to October time period. The model and observed results for ecosystem respiration (RESP) show that the peak value occurred during the middle of the summer with low values during the winter months. The RESP results generally followed the seasonal pattern of air temperature with maximum RESP and air temperatures (TEMP) during the middle of the summer, increasing RESP and TEMP during the spring, decreasing values during the fall and low values during winter months (temperature data not shown). Both the model results and observed 2005 NEE show that peak carbon uptake occurred in the spring (May), net carbon losses during the winter, and a short period of net carbon uptake in July and August. An inter-annual comparison of NEE for the variable and fix Amax versions of Daycent show the differences season responses of each model, with Fig. 6 show three consecutive years of simulations.

Comparison of cumulative annual average NEE show an average net uptake of carbon of 36.0 gC·m^‒2 for the observed data and 40.0 gC·m^‒2 for the model results. The results show that R² of the model vs observed for NEE and AET are lower than R² for GPP and RESP (R²>0.88 for GPP and RESP and R² of 0.66 and 0.67 for NEE and AET).

Discussion

Our results justify incorporating a photosynthetic sub-model and variable Amax into the DayCent model. The resulting DayCent-Photo simulated NEE, GPP, and ecosystem respiration (RESP) much better than the Base and the single Amax versions even after penalizing the model for including extra parameters. The results suggest that seasonality of GPP can be simulated by varying the capacity of canopy level photosynthesis. The biology behind the seasonality of GPP is certainly more complex than the variable Amax function we used, but the varying Amax maybe adequate enough to capture some aspects of the variability. The variable photosynthetic capacity has four known factors: 1) lower carbon uptake at the start of the growing season due to slow recovery of stomatal conductance (^{Monson et al., 2005}); 2) reduced carbon uptake due to increased temperature at mid-summer (^{Huxman et al., 2003}); 3) differences in assimilation rates of young versus old needles as canopy turns over during the growing season (^{Mohren and van de Veen, 1995}; ^{Marshall et al., 2001}; ^{McGarvey et al., 2004}; ^{Weiskittel et al., 2006}); and 4) an increase in leaf area index (^{Wang and Barrett, 2003}) resulting from needle growth during the growing season. The first and second factors are likely related to stomatal response to water availability, and the third and fourth factors may be related to the distribution of leaves within the layers of the canopy. The photosynthesis sub-model of DayCent-Photo is simple, with one-layer canopy and lumped parameter photosynthesis lacking stomatal conductance. This simplicity may have resulted in the single Amax version to poorly simulate the seasonality of GPP, but the simplicity was adequately augmented by the addition of the variable Amax function.

The model simulated NEE and especially GPP well across six years (Fig. 3), and likely can simulate them well for the near future given that climate change may shift the the beginning of seasonal growth by a few days per decade. For sub-alpine forests, how the seasonality will shift remains uncertain (^{Richardson et al., 2013}), but the timing of tree flowering has advanced by 1 to 5 days per decade (^{Ziello et al., 2009}). In boreal forests, spring leaf out has advanced by 0.3‒1.1 days per decade in Norway (^{Linkosalo et al., 2009}), and by 6.3 days per decade elsewhere in Europe (^{Delbart et al., 2008}). The model is able to simulate the advancing spring leaf out if the advancement is controlled by temperature. Leaf out and start of spring photosynthesis is a function of temperature in DayCent-Photo. The variable Amax function cannot predict how climate change affects the phenology of photosynthesis, but the phenology likely will change gradually for the near future, gradual enough for DayCent-Photo to predict NEE reasonably well. Seasonality of NEE may be due to temperature acclimation of photosynthesis, with colder sites slower to acclimate (^{Gea-Izquierdo et al., 2010}), suggesting that the optimization and validation should be site specific.

This study focused on the photosynthesis sub-model of DayCent, and the sub-model represents only a small portion of the total ecosystem processes that are simulated by the DayCent model. Since there are many papers that describe the setup and running of the DayCent model (^{Paustian et al., 1992}; Parton and Rasmussen, 1994; ^{Parton et al., 1998}), we limited our variability and analysis to only the new parameters added in the sub-model. The model performance could be further improved by adjusting other parameters, such as those controlling autotrophic and heterotrophic respiration rates. The DayCent-Photo model slightly overestimated NEE (more carbon uptake than observed, Fig. 3), as a result of overestimating GPP during the summer period and underestimating RESP in the onset and end of growing season. Climate change can hasten the date the forest transitions from carbon sink to source in autumn, primarily from increased RESP (^{Piao et al., 2008}), and more accurate simulation of RESP may be critical in predicting the carbon dynamics of the forest in the future.

The errors in NEE reported here do not include errors in the data. The eddy covariance technique typically underestimates RESP compared chamber-based measurements due to the influence of effects such as drainage and advection (^{Massman and Lee, 2002}). ^{Speckman et al. (2015)} showed that even with high wind speed RESP derived from the nighttime eddy covariance measurements was underestimated in a nearby forest in Wyoming. The errors associated with eddy covariance technique are estimated to be 30‒100 gC·m^‒2·yr^‒1 (^{Baldocchi, 2008}), and while DayCent-Photo overestimated the annual sum by less than 10 gC·m^‒2.

The DayCent-Photo model has recently been used to successfully simulate seasonal patterns in AET, NEE, GPP and RESP for deciduous forest, conifer forest, and grassland (^{Savage et al., 2013; Chen et al., 2016}). The results from all of these papers suggest that Amax has a similar seasonal patterns for grasslands and forest systems with Amax values being highest in the May to early June time period and low during the spring and fall months. The results show that including variable Amax greatly improved the ability of the model to simulate seasonal patterns in ecosystem variables and that observed seasonal patterns in GPP and NEE from the AmeriFlux sites can be used to successfully parameterize the seasonal changes in Amax. The recent results show that the DayCent-Photo model can be used to simulate ecosystem dynamics for different biomes (grasslands, conifer forests, and deciduous forest). We are currently expanding the use of the DayCent-Photo model to simulate seasonal changes in AET, NEE, GPP, and RESP for all of the major biomes (savanna, tell grass prairie and other forest types) in the world using data sets from AmeriFlux sites.

Conclusions

For this study we modified the DayCent model to include a photosynthetic sub-model with seasonally variable photosynthetic capacity. Compared to the base DayCent and single Amax versions, the inclusion of the photosynthesis sub-model and variable Amax was strongly supported by observations, even after the extra parameters were penalized for extra degrees of freedom. These results suggest that the seasonality of GPP and NEE in this system arises from the changes in photosynthetic capacity. The biological process behind the seasonality likely involves the complex biology of stomatal conductance, drought response, and leaf turnover.

By modeling the seasonal variation in the NEE flux through model optimization procedure, we extracted information about the necessary model structure needed to simulate the seasonality of photosynthesis. Extrapolating model structure to a biological process equation will be the next significant step in understanding the effect of Amax for different ecosystems. Once this is in place we can gain a further understanding the on the role variable Amax plays on ecosystems under various levels of stress. DayCent-Photo has been successfully used to simulate seasonal changes in ecosystem variables for grassland and forest ecosystems using AmeriFlux data sets to calibrate and test the model and has the potential to be successfully used to simulate ecosystem dynamics for all of the major ecosystems in the world.

References

[1] Aber J D, Federer C A (1992). A generalized, lumped-parameter model of photosynthesis, evapotranspiration and net primary production in temperate and boreal forest ecosystems. Oecologia, 92(4): 463–474

[2] Baldocchi D D (2008). Breathing of the terrestrial biosphere: lessons learned from a global network of carbon dioxide flux measurement systems. Aust J Bot, 56(1): 1–26

[3] Baldocchi D D (2003). Assessing the eddy covariance technique for evaluating carbon dioxide exchange rates of ecosystems: past, present and future. Glob Change Biol, 9(4): 479–492

[4] Baldocchi D D, Hincks B B, Meyers T P (1988). Measuring biosphere-atmosphere exchanges of biologically related gases with micrometeorological methods. Ecology, 69(5): 1331–1340

[5] Blanken P D, Monson R K, Burns S P, Turnipseed A A (2010). Data and Information for the US-NR1 Niwot Ridge Subalpine Forest AmeriFlux Site (LTER NWT1). AmeriFlux Management Project. Lawrence Berkeley National Laboratory, California

[6] Bourdeau P F (1959). Seasonal variations of the photosynthetic efficiency of evergreen conifers. Ecology, 40(1): 63–67

[7] Braswell B H, Sacks W J, Linder E, Schimel D S (2005). Estimating diurnal to annual ecosystem parameters by synthesis of a carbon flux model with eddy covariance net ecosystem exchange observations. Glob Change Biol, 11(2): 335–355

[8] Chen M, Parton W J, Adair E C, Asao S, Hartman M D, Gao W (2016). Simulation of the effects of photodecay on long-term litter decay using DayCent. Ecosphere, 7(12): e01631

[9] Chen M, Parton W J, Del Grosso S J, Hartman M D, Day K A, Tucker C J, Derner J D, Knapp A K, Smith W K, Ojima D S, Gao W (2017). The signature of sea surface temperature anomalies on the dynamics of semiarid grassland productivity. Ecosphere, 8(12): e02069

[10] Del Grosso S J, Parton W J, Mosier A R, Hartman M D, Brenner J, Ojima D S, Schimel D S (2001). Simulated Interaction of Carbon Dynamics and Nitrogen Trace Gas Fluxes Using the DAYCENT Model. In: Shaffer M J, Ma L W, Hansen S, eds. Modeling Carbon and Nitrogen Dynamics for Soil Management. Boca Raton: CRC Press, 303–332

[11] Del Grosso S J, Parton W J, Mosier A R, Ojima D S, Kulmala A E, Phongpan S (2000). General model for N₂O and N₂ gas emissions from soils due to dentrification. Global Biogeochem Cycles, 14(4): 1045–1060

[12] Del Grosso S J, Parton W J, Stohlgren T J, Zheng D L, Bachelet D, Prince S, Hibbard K, Olson R (2008). Global potential net primary production predicted from vegetation class, precipitation, and temperature. Ecology, 89(8): 2117–2126

[13] Delbart N, Picard G, Le Toan T, Kergoat L, Quegan S, Woodward I, Dye D, Fedotova V (2008). Spring phenology in boreal Eurasia over a nearly century time scale. Glob Change Biol, 14(3): 603–614

[14] Drake J E, Raetz L M, Davis S C, Delucia E H (2010). Hydraulic limitation not declining nitrogen availability causes the age-related photosynthetic decline in loblolly pine (Pinus taeda L.). Plant Cell Environ, 33(10): 1756–1766

[15] Fisher R A (1932). Inverse probability and the use of likelihood. Math Proc Camb Philos Soc, 28(03): 257–261

[16] Frey S D, Lee J, Melillo J M, Six J (2013). The temperature response of soil microbial efficiency and its feedback to climate. Nat Clim Chang, 3(4): 395–398

[17] Gea-Izquierdo G, Mäkelä A, Margolis H, Bergeron Y, Black T A, Dunn A, Hadley J, Paw U K T, Falk M, Wharton S, Monson R, Hollinger D Y, Laurila T, Aurela M, McCaughey H, Bourque C, Vesala T, Berninger F (2010). Modeling acclimation of photosynthesis to temperature in evergreen conifer forests. New Phytol, 188(1): 175–186

[18] Granda E, Scoffoni C, Rubio-Casal A E, Sack L, Valladares F (2014). Leaf and stem physiological responses to summer and winter extremes of woody species across temperate ecosystems. Oikos, 123(11): 1281–1290

[19] Guan M, Jin Z, Wang Q, Li Y, Zuo W (2014). Response of photosynthesis traits of dominant plant species to different light regimes in the secondary forest in the area of Qiandao Lake, Zhejiang, China. China Journal of Applied Ecology, 25: 1615–1622

[20] Hartman M D, Baron J S, Ewing H A, Weathers K C (2014). Combined global change effects on ecosystem processes in nine U.S. topographically complex areas. Biogeochemistry, 119(1−3): 85–108

[21] Helms J A (1965). Diurnal and seasonal patterns of net assimilation in Douglas-Fir, Pseudotsuga Menziesii (Mirb). Franco, as Influenced by Environment. Ecology, 46(5): 698–708

[22] Hilborn R, Mangel M (1997). The ecological detective: confronting models with data. Monogr Popul Biol, 28: 315

[23] Hurtt G C, Armstrong R (1996). A pelagic ecosystem model calibrated with BATS data. Deep Sea Res Part II Top Stud Oceanogr, 43(2−3): 653–683

[24] Huxman T E, Turnipseed A A, Sparks J P, Harley P C, Monson R K (2003). Temperature as a control over ecosystem CO₂ fluxes in a high-elevation, subalpine forest. Oecologia, 134(4): 537–546

[25] Johnson J B, Omland K S (2004). Model selection in ecology and evolution. Trends Ecol Evol, 19(2): 101–108

[26] Kelly R H, Parton W J, Hartman M D, Stretch L K, Ojima D S, Schimel D S (2000). Intra-annual and interannual variability of ecosystem processes in shortgrass steppe. Journal of Geophysical Research: Atmospheres, 105(D15): 20093–20100

[27] Li Z, Li X, Rubert-Nason K F, Yang Q, Fu Q, Feng J, Shi S (2018). Photosynthetic acclimation of an evergreen broadleaved shrub (Ammopiptanthus mongolicus) to seasonal climate extremes on the Alxa Plateau, a cold desert ecosystem. Trees (Berl), 32(2): 603–614

[28] Linkosalo T, Häkkinen R, Terhivuo J, Tuomenvirta H, Hari P (2009). The time series of flowering and leaf bud burst of boreal trees (1846-2005) support the direct temperature observations of climatic warming. Agric Meteorol, 149(3–4): 453–461

[29] Luyssaert S, Ciais P, Piao S L, Schulze E D, Jung M, Zaehle S, Schelhaas M J, Reichstein M, Churkina G, Papale D, Abril G, Beer C, Grace J, Loustau D, Matteucci G, Magnani F, Nabuurs G J, Verbeeck H, Sulkava M, van der WERF G R, Janssens I A (2010). The European carbon balance. Part 3: forests. Glob Change Biol, 16(5): 1429–1450

[30] Luyssaert S, Schulze E D, Börner A, Knohl A, Hessenmöller D, Law B E, Ciais P, Grace J (2008). Old-growth forests as global carbon sinks. Nature, 455(7210): 213–215

[31] Marshall J D, Rehfeldt G E, Monserud R A (2001). Family differences in height growth and photosynthetic traits in three conifers. Tree Physiol, 21(11): 727–734

[32] Martinez K A, Fridley J D (2018). Acclimation of leaf traits in seasonal light environments: Are non-native species more plastic? J Ecol, 20: 207–216

[33] Massman W J, Lee X (2002). Eddy covariance flux corrections and uncertainties in long-term studies of carbon and energy exchanges. Agric Meteorol, 113(1−4): 121–144

[34] McGarvey R C, Martin T A, White T L (2004). Integrating within-crown variation in net photosynthesis in loblolly and slash pine families. Tree Physiol, 24(11): 1209–1220

[35] Mohren G M J, van de Veen J R (1995). Forest growth in relation to site conditions. Application of the model forgro to the Solling spruce site. Ecol Modell, 83(1−2): 173–183

[36] Monson R K, Sparks J P, Rosenstiel T N, Scott-Denton L E, Huxman T E, Harley P C, Turnipseed A A, Burns S P, Backlund B, Hu J (2005). Climatic influences on net ecosystem CO₂ exchange during the transition from wintertime carbon source to springtime carbon sink in a high-elevation, subalpine forest. Oecologia, 146(1): 130–147

[37] Monson R K, Turnipseed A A, Sparks J P, Harley P C, Scott-Denton L E, Sparks K, Huxman T E (2002). Carbon sequestration in a high-elevation, subalpine forest. Glob Change Biol, 8(5): 459–478

[38] Moore D J P, Hu J, Sacks W J, Schimel D S, Monson R K (2008). Estimating transpiration and the sensitivity of carbon uptake to water availability in a subalpine forest using a simple ecosystem process model informed by measured net CO₂ and H₂O fluxes. Agric Meteorol, 148(10): 1467–1477

[39] Papale D, Valentini R (2003). A new assessment of European forests carbon exchanges by eddy fluxes and artificial neural network spatialization. Glob Change Biol, 9(4): 525–535

[40] Parton W J, Hanson P J, Swanston C, Torn M, Trumbore S E, Riley W, Kelly R (2010). ForCent model development and testing using the enriched background isotope study experiment. J Geophys Res, 115(G4): G04001

[41] Parton W J, Hartman M, Ojima D, Schimel D (1998). DAYCENT and its land surface submodel: description and testing. Global Planet Change, 19(1−4): 35–48

[42] Parton W J, Rasmussen P E (1994). Long-term effects of crop management in wheat/fallow: II. CENTURY model simulations. Soil Sci Soc Am J, 58(2): 530–536

[43] Parton W, Holland E A, Del Grosso S J, Hartman D, Martin M, Mosier A, Ojima D S, Schimel D S (2001). Generalized model for NO_x and N₂O emissions from soils. J Geophys Res, 106(D15): 17403–17419

[44] Paustian K, Parton W J, Persson J (1992). Modeling soil organic matter in organic-amended and nitrogen-fertilized long-term plots. Soil Sci Soc Am J, 56(2): 476–488

[45] Piao S, Ciais P, Friedlingstein P, Peylin P, Reichstein M, Luyssaert S, Margolis H, Fang J, Barr A, Chen A, Grelle A, Hollinger D Y, Laurila T, Lindroth A, Richardson A D, Vesala T (2008). Net carbon dioxide losses of northern ecosystems in response to autumn warming. Nature, 451(7174): 49–52

[46] Rastetter E B, Aber J D, Peters D P C, Ojima D S, Burke I C (2003). Using mechanistic models to scale ecological processes across space and time. Bioscience, 53(1): 68

[47] Reichstein M, Falge E, Baldocchi D, Papale D, Aubinet M, Berbigier P, Bernhofer C, Buchmann N, Gilmanov T, Granier A, Grunwald T, Havrankova K, Ilvesniemi H, Janous D, Knohl A, Laurila T, Lohila A, Loustau D, Matteucci G, Meyers T, Miglietta F, Ourcival J M, Pumpanen J, Rambal S, Rotenberg E, Sanz M, Tenhunen J, Seufert G, Vaccari F, Vesala T, Yakir D, Valentini R (2005). On the separation of net ecosystem exchange into assimilation and ecosystem respiration: review and improved algorithm. Glob Change Biol, 11(9): 1424–1439

[48] Richardson A D, Keenan T F, Migliavacca M, Ryu Y, Sonnentag O, Toomey M (2013). Climate change, phenology, and phenological control of vegetation feedbacks to the climate system. Agric Meteorol, 169: 156–173

[49] Ryan M G, Waring R H (1992). Maintenance respiration and stand development in a young subalpine lodgepole pine forest. Ecology, 73: 2100–2108

[50] Sacks W J, Schimel D S, Monson R K (2007). Coupling between carbon cycling and climate in a high-elevation, subalpine forest: a model-data fusion analysis. Oecologia, 151(1): 54–68

[51] Sacks W J, Schimel D S, Monson R K, Braswell B H (2006). Model-data synthesis of diurnal and seasonal CO₂ fluxes at Niwot Ridge, Colorado. Glob Change Biol, 12(2): 240–259

[52] Savage K E, Parton W J, Davidson E A, Trumbore S E, Frey S D (2013). Long-term changes in forest carbon under temperature and nitrogen amendments in a temperate northern hardwood forest. Glob Change Biol, 19(8): 2389–2400

[53] Schimel D (1995). Terrestrial ecosystems and the carbon cycle. Glob Change Biol, 1(1): 77–91

[54] Speckman H N, Frank J M, Bradford J B, Miles B L, Massman W J, Parton W J, Ryan M G (2015). Forest ecosystem respiration estimated from eddy covariance and chamber measurements under high turbulence and substantial tree mortality from bark beetles. Glob Change Biol, 21(2): 708–721

[55] Tang X, Wang X, Wang Z, Liu D, Jia M, Dong Z, Xie J, Ding Z, Wang H, Liu X (2013). Influence of vegetation phenology on modelling carbon fluxes in temperate deciduous forest by exclusive use of MODIS time-series data. Int J Remote Sens, 34(23): 8373–8392

[56] Tang X, Wang Z, Liu D, Song K, Jia M, Dong Z, Munger J W, Hollinger D Y, Bolstad P V, Goldstein A H, Desai A R, Dragoni D, Liu X (2012). Estimating the net ecosystem exchange for the major forests in the northern United States by integrating MODIS and AmeriFlux data. Agric Meteorol, 156: 75–84

[57] Turnipseed A A, Anderson D E, Blanken P D, Baugh W M, Monson R K (2003). Airflows and turbulent flux measurements in mountainous terrain. Part 1. Canopy and local effects. Agric Meteorol, 119(1−2): 1–21

[58] Turnipseed A A, Anderson D E, Burns S, Blanken P D, Monson R K (2004). Airflows and turbulent flux measurements in mountainous terrain: Part 2: Mesoscale effects. Agric Meteorol, 125(3−4): 187–205

[59] Turnipseed A A, Blanken P D, Anderson D E, Monson R K (2002). Energy budget above a high-elevation subalpine forest in complex topography. Agric Meteorol, 110(3): 177–201

[60] Urban O, Holub P, Klem K (2017). Seasonal courses of photosynthetic parameters in sun- and shade-acclimated spruce shoots. Beskydy, 10(1−2): 49–56

[61] Wang Y P, Baldocchi D, Leuning R, Falge E, Vesala T (2007). Estimating parameters in a land-surface model by applying nonlinear inversion to eddy covariance flux measurements from eight FLUXNET sites. Glob Change Biol, 13(3): 652–670

[62] Wang Y P, Barrett D J (2003). Estimating regional terrestrial carbon fluxes for the Australian continent using a multiple-constraint approach I. Using remotely sensed data and ecological observations of net primary production. Tellus B Chem Phys Meterol, 55: 270–289

[63] Weiskittel A R, Maguire D, Garber S M, Kanaskie A (2006). Influence of Swiss needle cast on foliage age-class structure and vertical foliage distribution in Douglas-fir plantations in north coastal Oregon. Can J Res, 36(6): 1497–1508

[64] Zhang Y J, Holbrook N M, Cao K F (2014). Seasonal dynamics in photosynthesis of woody plants at the northern limit of Asian tropics: potential role of fog in maintaining tropical rainforests and agriculture in Southwest China. Tree Physiol, 34(10): 1069–1078

[65] Zhang Y J, Sack L, Cao K F, Wei X M, Li N (2017). Speed versus endurance tradeoff in plants: leaves with higher photosynthetic rates show stronger seasonal declines. Sci Rep, 7(1): 42085

[66] Ziello C, Estrella N, Kostova M, Koch E, Menzel A (2009). Influence of altitude on phenology of selected plant species in the Alpine region (1971-2000). Clim Res, 39: 227–234

[67] Zobitz J M, Moore D J P, Sacks W J, Monson R K, Bowling D R, Schimel D S (2008). Integration of process-based soil respiration models with whole-ecosystem CO₂ measurements. Ecosystems (N Y), 11(2): 250–269

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag GmbH Germany, part of Springer Nature