Representing Grasslands Using Dynamic Prognostic Phenology Based on Biological Growth Stages: 1. Implementation in the Simple Biosphere Model (SiB4)

Grasslands grow in a sequence of seasonal growth stages that respond to both climate and weather, and these relationships can be used to establish a strategy for predicting plant phenology. Current plant states (phenophase) can be represented as one of established growth stages that dictate carbon allocation and leaf photosynthetic capacity. Calculating daily phenophases from climate and environmental relationships allows for sequential growth stages (i.e., well‐defined seasonal cycles with a single growth period) or dynamic growth stages (i.e., multiple growth periods during a growing season). Senescence results from biomass mortality in response to environmental conditions. This approach uses a single mechanistic framework to represent grassland ecology, removing the dependence on satellite‐based vegetation indices and individual site tuning of parameters. Rather than being specified, a variety of properties emerge, from allometric relationships such as root‐shoot ratios, to behavior across moisture gradients, to interannual variability in growing season lengths, carbon stores, and land surface fluxes. Using dynamic phenology stages to link biophysical and biogeochemical processes provides a mechanism to predict self‐consistent land‐atmosphere exchanges of carbon, water, energy, radiation, and momentum, as well as carbon storage in cascading pools of biomass; and describing these processes in a mathematically determinate model makes them clear, testable, and usable for predictions. This paper describes this new phenology method as it is implemented in the Simple Biosphere Model Version 4 (SiB4), and a companion paper evaluates this method at grassland sites worldwide.


Introduction
Phenology in terrestrial ecosystems, which refers to periodic vegetation events that are driven by environmental conditions, plays an essential role in biogeochemical processes, from determining local land surface properties to impacting regional land-atmosphere interactions, driving the terrestrial carbon cycle, and influencing global climate change (Dahlin et al., 2015;Friedl et al., 2014;Kim et al., 2015;Schwartz, 2013). Plant phenology, which incorporates visible seasonal events while integrating underlying biological functionality, can be easily seen and understood, making it a powerful indicator of climate change; however, since phenology provides a holistic view of biological mechanisms, it is also a powerful tool for scientists.
Since there is a critical link between phenology and terrestrial biochemical, biophysical, and spectral properties, studying phenology not only leads to further understanding of ecosystem processes but also can be used to monitor and understand biospheric responses to global climate change (Fu et al., 2014;Gu et al., 2003;Richardson et al., 2013).
Vegetation phenology is a critical component of terrestrial biosphere models (TBMs), which are essential tools for a variety of applications, from predicting land-atmosphere exchanges to improving our understanding of terrestrial ecosystems and predicting impacts of climate change. Although it is well known that phenology is specific to individual species, TBMs provide an opportunity to explore ecosystem-level seasonal patterns and capture broad vegetational responses (Gu et al., 2003). Accurately simulating phenology is essential in TBMs because seasonality controls numerous aspects of plant growth and development, which then regulates ecosystem processes and land-atmosphere fluxes, directly impacting terrestrial biogeochemical and hydrological cycles (Eamus et al., 2016;Friedl et al., 2014;Pau et al., 2011).
To simulate plant phenology, TBMs employ two different methods: using satellite observations directly or predicting phenology using process-based models. While prescribing the phenology from satellite observations takes advantage of global remotely sensed data, these models are limited to simulating time periods when satellite observations are available (Medvigy et al., 2009;. In contrast, TBMs that model phenology have the advantage of predicting the vegetation state at the expense of introducing empirical or calibrated parameters (Chuine et al., 2013;Kim et al., 2015;Lawrence et al., 2011).
Due to the strong applicability to climate change research, interest in including process-based phenology models in TBMs has been expanding. One common approach to predict phenology used by TBMs is a rule-based strategy to predict phenological transitions. Environmental conditions are prominent causes of plant phenology, and of these temperature is strongly related to leaf-out and growing season onset (Schwartz, 2013). In using temperature to time phenological transitions, one technique is to use a growing-degree day approach White et al., 1997). For simulating crops, the growing-degree day determines the developmental stage, which is used to adjust carbon allocation and vegetation characteristics as the growth progresses throughout the season (e.g., Lokupitiya et al., 2009). For simulating natural vegetation, global TBMs use cumulative growing-degree days to identify leaf onset (Lawrence et al., 2011;Krinner et al., 2005). To improve phenological predictions, rule-based models employ chilling requirements for leaf onset (e.g., Kaduk & Los, 2011;Picard et al., 2005), add day of year as an additional parameter (e.g., Fisher et al., 2007), or incorporate soil water and photoperiod for the timing of senescence (e.g., Borchert et al., 2005;Dahlin et al., 2015Dahlin et al., , 2017Lawrence et al., 2011).
An alternative to predicting phenological transitions is to continuously characterize leaf phenology. To do this, Jolly et al. (2005) introduced a combined growing season index that utilizes thresholds of temperature, light, and humidity to predict the vegetation state. Stöckli et al. (2008Stöckli et al. ( , 2011 used this approach combined with data assimilation to fit parameters at every grid cell and predict present-day phenology. Forkel et al. (2014) adopted this approach in a dynamic vegetation model, replacing the humidity requirement with water availability and adding a heat stress limiting function. In a different approach that still continuously characterizes leaf phenology, Caldararu et al. (2014) developed a mechanistic model of leaf phenology using the hypothesis that phenology is a strategy for optimal carbon gain at the canopy level, making leaf adjustments in response to light, temperature, and soil moisture.
Although work is ongoing to improve phenology predictions, TBMs have systematic errors in the simulation of phenology (Migliavacca et al., 2012). A model intercomparison study evaluated the representation of phenology and the associated seasonality of ecosystem-scale CO 2 exchange , finding that existing models are unable to accurately simulate the timing of both the beginning and ending of the growing season and highlighting the need for improved understanding and representation of environmental controls on vegetation phenology. Additionally, a study by Keenan et al. (2012) found that models underrepresent the variability in spring leaf-out, leading to systematic errors in carbon fluxes. When used in a global climate model, TBMs with prognostic vegetation phenology add the capability to represent interannual variability; however, mean biogeophysical simulations tend to be degraded due to the added degree of freedom (Lawrence et al., 2011(Lawrence et al., , 2012. One of the difficulties in predicting phenology in TBMs comes from the fact that the ecological processes underlying observed phenological patterns at ecosystem scales is still uncertain (Chuine et al., 2013; 10.1029/2018MS001540 Liang & He, 2017;Melaas et al., 2016;Pau et al., 2011;Richardson et al., 2013). Satellite remote sensing provides broad spatial and temporal coverage of terrestrial vegetation; however, the attribution of observed land surface phenology to biophysical processes remains lacking (Liang & He, 2017;Richardson et al., 2017). Increasing observational data and real-time monitoring combined with cross-scale data integration and modeling provides an opportunity for improving the understanding of the relationships between phenology patterns and ecological processes from local to regional to global scales (Richardson et al., 2017;Morellato et al., 2016).
We developed a prognostic phenology strategy that uses a combination of predicting rule-based transitions while continuously calculating leaf phenology and the vegetation state in an innovative hybrid approach. This method continuously predicts the current phenological stage (phenophase) by creating a phenology index (PI) that is determined from three elements: (1) a growth index (PS GI ) that estimates the benefit of growing new leaves; (2) a weather potential (PS WX ) that calculates vegetation stress based on current temperature and water availability; and (3) a day length index (PS DL ) to obtain developmental responses to photoperiod. The phenophase is used to determine the carbon allocation and maximum photosynthesis rate. Sencescence is simulated from both continual turnover (with rates that respond to assimilation rate and temperature) and increasing transfer fractions (with fractions that depend on daylength, temperature, and water availability). The combination of these elements results in the simulation of key phenological events (such as leaf-out and senescence) while introducing an innovative, dynamic approach. By unifying phenology, carbon pools, and land-atmosphere exchanges, the Simple Biosphere Model (SiB4) predicts self-consistent land surface properties, vegetation characteristics; soil moisture content; carbon storage in cascading pools of biomass; and land-atmosphere exchanges of carbon, water, energy, radiation, and momentum. This paper describes the new phenology strategy and its implementation in SiB4 Sellers et al., 1986;Sellers, Los, et al., 1996;. In representing grasslands, SiB4 uses three plant functional types (PFTs): arctic C3 (C3A), non-arctic C3 (C3G), and C4 (C4G). All parameters used in this methodology are PFT-specific. Section 2 describes the use of phenological stages to predict leaf-out and vegetation states throughout the growing season, and section 3 describes the use of turnover and litterfall processes to predict senescence. Since grazing is important in grassland ecosystems, section 4 briefly describes the implementation of grazing, and section 5 describes how these processes are combined to fully predict the terrestrial carbon cycle. While a thorough evaluation of this approach is in a companion paper, section 6 provides an example of the terrestrial carbon cycle at a temperate grassland site. Finally, section 7 discusses the utility of this approach and provides suggestions for future applications. Supporting information provides justification for the use of 10-day running means (section S1), grazing equations (section S2), and a step-by-step description and illustration of the calculation sequence in SiB4 (section S3), along with tables of the new variables, parameter values, and figures demonstrating the relationships used in the phenology model.
Since this paper only shows results at one site merely to illustrate the methodology, a companion paper (Haynes et al., 2019) thoroughly investigates the terrestrial carbon cycle of 55 grassland sites worldwide, using what we believe to be an unprecidented collection of grassland data. The companion paper uses remotely sensed leaf area index (LAI), eddy covariance carbon and energy fluxes at towers, and annual net primary productivity at long-term in situ sites to evaluate diurnal, seasonal, annual, and interannual variability predicted by SiB4 compared to these data sets while simultaneously investigating large-scale grassland relationships spanning climate zones, including the impact of precipitation gradients and drought on phenology, aboveground vegetation, and productivity. Haynes et al. (2019) shows that SiB4 predicts LAI, carbon, and energy fluxes with root-mean-square errors <15% and individual biases <10% at all available sites. SiB4 is able to predict observed large-scale grassland behavior across precipitation gradients and in response to drought with nearly identical relationships to those seen in the observations.

Dynamic Prognostic Phenology
Plant development responds to seasonal and interannual climatic variations (Jarvis, 1976;Larcher, 2002;Stöckli et al., 2008). Seasonal development occurs via physiological processes that cause visual changes in biomass. All plants pass through the main phenological events in an annual cycle, and these visual characteristics can be categorized into phenological stages. While phenology progressions are rhythmical, they are never identical nor exactly periodic; instead, the rates at which phenological changes occur depend on environmental conditions and differ not only between vegetation types but also in time (Schwartz, 2013). Phenological events are results of internal plant processes, thus continuously predicting the phenophase provides a mechanism to procure not only the visible changes in vegetation, but also the observed changes in physiology and biogeochemistry.
In order to respond to climate and weather conditions, plant growth must depend on environmentally responsive mechanisms from a combination of contributions, including air temperature, precipitation, humidity, radiation, and soil conditions (Chen, 2003). Phenophases can be predicted from climatological and environmental relationships. For natural vegetation, determining the phenophase is challenging as there are numerous confounding issues, including regions with no clear growing season and plants that have multiple growing seasons (Reed et al., 2003).

Phenological Stages
To predict grassland phenology, SiB4 uses dynamic phenophases that respond to the leaf state and environmental conditions. The overarching idea is that vegetation progresses through phenology stages during a growing season; however, the phenophase must allow plants to progress through the stages at widely different rates or to return to growth stages after midseason stress. For example, desert environments with short rainy seasons rapidly progress through the stages, while tropical grasslands are never dormant and rarely experience full leaf-out since they maintain leaves year-round.
To represent grasslands as they develop and mature through the growing season, this dynamic phenology scheme defines five seasonal stages (nstage = 5): 1. Leaf-Out: Leaf growth begins from stored carbon, since the canopy cannot support photosynthesis. This stage begins when the environmental conditions become suitable for photosynthesis, which is determined from daylength, soil moisture, and temperature. 2. Growth: Leaf growth is promoted from carbon allocation, since the canopy benefits by intercepting more light. This stage begins when the canopy is large enough to support photosynthesis. 3. Maturity: Leaf amount is sustained and allocation is balanced between canopy maintenance, root growth, and product production (seeds, fruit, and flowers). This stage begins when the plant no longer benefits from growing new leaves. 4. Stress: All carbon is allocated to the roots for storage and leaf photosynthetic capacity (V Max ) is reduced.
This stage begins when there is a penalty for growth (negative cost for new leaves). 5. Dormancy: Canopy is inert, there is no photosynthesis, and autotrophic respiration continues at a low rate. This stage begins when there are no longer leaves in the canopy, and continues until phenological triggers initiate a transition into leaf-out.
The remote sensing community supports the use of five stages for natural vegetation, as time series of satellite data are used to identify key seasonal transitions in vegetation activity between four stages: greenup, maturity, senescence, and dormancy (Hanes et al., 2013;Zhang et al., 2003). After performing sensitivity studies, we separated greenup into two stages (leaf-out and growth) because we found that five stages optimally captures the leaf-out and senescence, while providing sufficient options to simulate carbon flux and pool day-to-day, seasonal, and interannual variability using a mechanistic approach.

Phenology Stage Selection
The timing and length of each phenology stage is not prescribed, but rather is diagnosed from assimilation rate, climate, day length, leaf pool size, and plant stress from the associated weather and environmental conditions. In order to quantify the shifting between phenology stages, this scheme uses a phenology index (PI) such that where PS GI is a growth index, PS Wx is a weather potential, and PS DayL is a day length potential. PI will be high when it is beneficial to grow new leaves and when the environmental conditions support photosynthesis. In contrast, PI will be 0 either when there is no benefit in growing new leaves or when conditions are not suitable for photosynthesis. To distinguish between the phenology stages, SiB4 applies prescribed threshold values (PI Thresh ) to PI.

Growth Index (PS GI )
The growth index determines the benefit of growing new leaves as opposed to allocating carbon belowground. We developed this index using the light-use efficiency curve, radiation saturation effect, and climatic sustainability. When the vegetation coverage is low, the addition of each leaf has a significant impact on the photosynthesis rate; however, as the canopy closes and approaches the maximum aboveground biomass that can be supported by the climate, adding more leaves is no longer beneficial. To mimic this behavior, the growth index changes as necessary, starting at its maximum value when the vegetation is minimal (highest potential for growth) and decreasing to its minimal value with a full canopy.
The growth index (PS GI ) linearly decreases with increasing LAI such that where ClimLAI Min and ClimLAI Max are the diagnosed minimum and maximum LAI for a given climate and PSG Max and PSG Min are parameters.
In order to determine the canopy fullness supported by climatic conditions, this method uses climatological potentials that are 10-year running means initialized during spin-up. ClimLAI Min and ClimLAI Max are calculated such that where ClimP is a climatological suitability index and CL C , CL L , and CL G are parameters.
ClimP can exponentially or linearly scale climatic water availability to cover environmental conditions specific to each vegetation type and is given by where ClimW is a climatological water availability index and CP A , CP B , CP C , CP D , CP Min , and CP Max are parameters. ClimP is minimal in climates that cannot support vegetation and increases with increasing water availability.
ClimW provides a standardized metric to determine the climatological amount of soil water available to plants, since water availability is a key determinant of plant productivity. For grasslands, ClimW is the 10-year running mean of the root-weighted fractional plant available water (PAW ). PAW is defined as where RootF i is the root fraction in soil layer i, VL i is the volume of liquid water in soil layer i, FC is the field capacity, and WP is the wilting point. Since PAW only includes liquid water, this metric is lowered by times and/or layers of soil that are frozen. By averaging PAW over 10 years, ClimW provides a fractional metric for the climatological liquid water obtainable to plants through their roots, with unity representing unrestricted access to water and 0 indicating no water is available.
Relationships between ClimP and ClimW are shown in Figure S1. C4 grasslands use the exponential functionality (using CP A and CP B ) to rapidly increase ClimP for water availability values in the midrange of ClimW. For these grasslands, the scale of ClimP stretches to values between 0.1 to 2, since small differences in water availability can produce large differences in vegetation greenness and productivity. C3 grasslands use the linear relationship between ClimP and ClimW. Since temperate C3 grasslands cover broader climate conditions, they span the entire range of ClimW and expand the range from one to two. In contrast, arctic grasslands are honed in on the lower ClimW values, since much of the year the ground is frozen.
Using PFT-specific relationships, this method uses site-specific long-term water availability to create climatological minimum and maximum bounds of LAI that can then be compared to current LAI values, determining PS GI . Example PS GI values are shown in Figure S2, where each colored line represents different site-specific ClimP values determining the LAI bounds. At the beginning of the growing season when the LAI is minimal, PS GI is at its maximum potential (PS GI,Max ). As the LAI increases, the benefit of adding more leaves decreases due to increasing respiratory costs or decreasing resources, such as limited light or water. PS GI linearly decreases with increasing LAI until it reaches its minimum value (PS GI,Min ) when LAI surpasses ClimLAI Max . As PS GI decreases, PI decreases, advancing the phenology stage.

Weather Potential (PS Wx )
The phenology weather potential, PS Wx , utilizes real-time weather conditions and soil moisture to represent the potential photosynthesis. PS Wx is unity at the start of the growing season when environmental conditions are suitable for photosynthesis, ensuring rapid growth immediately at the onset of growth. As the growing season progresses and environmental conditions required for photosynthesis deteriorate, PS Wx decreases, until it reaches zero under highly stressed conditions when photosynthesis is inhibited.
The equation used to calculate PS Wx is where F Wx,RM is combined potential that is sensitive to soil moisture, humidity, and temperature and F Wx,SM is the seasonal maximum of F Wx,RM . F Wx,RM is given by where F E,RM and F WA,RM are 10-day running mean potentials representing environmental conditions and water availability. Since the overarching purpose of the phenology scheme is to follow the natural plant progression through the growing season, 10-day running means are used to avoid spontaneous reactions in plant growth due to short-lived weather systems while still responding to synoptic and seasonal-scale anomalies. The running mean length is an input parameter and can be changed as the user sees appropriate; however, SiB4 uses a 10-day running mean for all grasslands, and sensitivity studies justifying this choice are provided in section S1.
F E is a combined potential from three environmental conditions, such that All three conditions are potentials that equal unity under optimal conditions and 0 in adverse environmental conditions (Sellers, 1987). The F LH potential is a function of water vapor pressure and represents leaf humidity stress (Sellers et al., 1992). The F RZ potential represents root-zone water stress and can be summarized as follows: if the total column soil moisture is above the field capacity, then there is maximum potential and no stress; in contrast, if the soil moisture is below the wilting point, then the plant is totally stressed and thus has zero potential . The F T factor is a function of temperature and incorporates high-and low-temperature stress Jarvis, 1976).
To intensify the effects of water availability, F WA,RM is the 10-day running mean of the root-weighted fraction of PAW in the top three soil layers (PAW Top ), which is given by where RootF 3 is the total rooting fraction in the top three soil layers.

Day Length Potential (PS DayL )
Photoperiod is a dominant cue to initiate senescence processes in the arctic (Ernakovich et al., 2014;Estiarte & Peñuelas, 2015;Pau et al., 2011). To progress the phenology as days get shorter, as well as to prohibit leaf-out during mild weather events during the winter, this strategy uses a photoperiod potential (PS DayL ).
As the days shorten, PS DayL decreases, lowering PI and advancing the phenology stage. The equations for PS DayL are  Figure S4.

Growing Season Start
Unlike crops that have a definitive growing season that begins with a planting date, the beginning of the growing season for grasslands is not as straightforward and depends on climate: temperate grasslands may have strong seasonality with a dormancy period during the winter, while tropical grasslands may continue growing all year. Adding even more complexity, temperate grasslands may have two growing periods within a single year, with a brief period of browning in-between when the environmental conditions are too hot and dry to support growth. To mimic natural behavior, the phenology scheme needs to be dynamic, with the capacity to return to the leaf-out or growing stages when necessary, while not triggering leaf-out midseason for vegetation with a clearly defined single seasonal cycle upon experiencing midseason stress.
To accomplish this, one of the phenological potentials, PS Wx can be reset in order to help obtain the leaf-out stage when warranted. With PS Wx being a ratio, it is reset to unity by setting the seasonal maximum value F E,SM equal to the current value F E,RM . For seasonal vegetation, this ensures rapid growth from a minimal leaf pool when day lengths are not substantially decreasing. However, for nonseasonal vegetation that may not ever obtain the conditions required for resetting this potential, leaf-out or growth can still be obtained if the conditions support rapid growth.
To determine when F E,SM is reset, SiB4 uses two flags. The first is the growing season suitability flag (FlagG), which is true when conditions are suitable for photosynthesis. Following the growing season start drivers used by Jolly et al. (2005), SiB4 uses temperature, soil moisture, and day length requirements to begin a growing season. For temperature, the daily mean temperature (TM D ) must be greater than a specified temperature (parameter FGT Min ) for a specific number of days (parameter FGT Len ). For soil moisture, SiB4 uses the root-weighted total available fraction of water in the top three soil layers (TAW Top ), which is given by where VI i is the volume of ice in soil layer i. Whereas PAW is penalized during time periods of frozen soil, TAW includes both liquid and ice water contributions and thus is a metric for total water availability. The top three soil layers (26 mm depth) contain ∼75% of simulated grassland roots; thus, these layers can be used to indicate water availability while responding quickly to precipitation events given their relatively shallow depth. To start the growing season, TAW Top must be above a specified value (parameter FGW Min ) for a specific number of days (parameter FGW Len ). For light, the day length must be longer than a specified minimum length. When the days are getting longer, the required day length is given by the parameter FGD MinI .

10.1029/2018MS001540
When the days are getting shorter, the required day length is given by the parameter FGD MinD . If all three of these conditions are met, then FlagG is true, otherwise FlagG is false.
The second flag used to determine the start of the growing season is an assimilation flag (FlagA), which is set to true when the assimilation rate becomes minimal at the end of the growing season. The assimilation flag makes use of an assimilation potential (F A ) and is true when where FlagA Thresh is a parameter, AssimD RM is a 10-day running mean of the daily mean assimilation rate AssimD, and AssimD SM is the seasonal maximum of AssimD RM . AssimD SM is reset to AssimD RM at the beginning of the growing season, when F E,SM is also reset. The occurrence of both FlagA and FlagG indicates that the vegetation can commence leaf-out to start a new growing season.

Phenology and Physiology Interactions
We hypothesize that two physiological processes are directly dependent on phenophase: carbon allocation and leaf photosynthetic capacity. Carbon allocation is a critical physiological process that determines the fate of assimilated carbon and affects plant growth (Xia et al., 2017). A compilation of research suggests that allocation to the leaf pool decreases with advancing phenology stage, with a shift in allocation from aboveground to belowground through the growing season ( . This progression makes sense intuitively as well: at the beginning of the growing season when the leaf pool is small, plants are limited by the light they collect and thus allocate resources to growing leaves. At the end of the season plants allocate more resources belowground to the roots for two reasons: (1) in many ecosystems water becomes a limiting resource (e.g., Blair et al., 2014;Rishmawi et al., 2016;Sala et al., 1988) and (2) roots store carbon required for growth at the beginning of the next growing season.
Using the philosophy that plants adjust their allocation to improve their most limiting resource, SiB4 includes seasonally varying allocation by including stage-specific allocation fractions (parameter AllocP nlpool, nstage ) that specify the fraction of photosynthesized carbon allocated to each live pool for each phenology stage. Allocation to the leaf pool in the case of light limitation is captured in the phenology stage selection: if the environmental conditions support growth and if adding leaves would provide an improvement to photosynthesis at a reasonable cost, then the phenology stage shifts to the leaf-out or growth stages with high leaf allocation fractions. Moisture sensitivity is also included in the phenology stage determination: if the vegetation is moisture stressed, then the phenophase will shift toward later stages that have more allocation to the roots. While the phenophase captures allocation changes throughout the season, complimentary weather-based adjustments provide a way to capture day-to-day variability, with adjustments being made based on water and temperature conditions (Friedlingstein et al., 1999;Schaefer et al., 2008).
Second, plants can change their leaf photosynthetic capacity (maximum carboxylation rate of rubisco, V Max ) throughout the season. V Max is an important control on the photosynthetic rate, changing the amount of carbon available for growth and maintenance and having significant impacts on the carbon cycle, energy exchanges, and hydrology (Baldocchi et al., 2001;Bauerle et al., 2012;Groenendijk et al., 2011;Wilson & Baldocchi, 2000). Increasing evidence indicates that large seasonal variations in V Max occur in terrestrial ecosystems (Bauerle et al., 2012;Groenendijk et al., 2011;Houborg et al., 2013;Zhou et al., 2014); and specifying V Max by phenophase ties it to seasonal development and provides a mechanistic way to change it throughout the growing season. V Max increases through the growth stage, then decreases as the leaves age (Restrepo-Coupe et al., 2013;Saleska et al., 2014;Suzuki et al., 2010;Wu et al., 2016;Xu et al., 2017).
Implementing prognostic phenology using dynamic stages provides a powerful mechanism to integrate ecosystem processes, linking carbon pools (storage), land-atmosphere (energy) fluxes, and biosphere-atmosphere (carbon) exchanges. Phenophases can be determined from environmental conditions and used to connect physiology and biogeochemistry, and Figure 1 provides an illustration of this integration. Using phenophase to change plant physiology shifts the focus from prescribing physical properties, which are highly variable in space and time, to predicting them. Rather than using empirical allometric relationships, ecosystem characteristics emerge from processes that vary with weather and climate. For example, rather than using a predefined root-to-shoot ratio, this is an emergent property in this model, which is bounded by climate and weather interactions rather than by specified parameters.

Turnover
Similar to SiB-CASA (Schaefer et al., 2008), pool turnover is tied to respiration; however, SiB4 introduces new scaling coefficients that are described in this section. To capture the loss of carbon from any live pool (lp) that is transferred to a dead carbon pool, SiB4 uses the equation where LossT lp,t is the loss due to pool turnover at time step t, E lp is the pool respiration efficiency, C lp is the pool carbon, lp is the pool turnover time, and dt is the time step length. M A,t , M F,t , and M H,t are the assimilation rate, freeze inhibition, and high temperature scaling coefficients at time step t, respectively. The base turnover time is the longest for the product pool based on the premise that seeds and flowers are carbon intensive to produce yet require less energy to maintain, and it is the shortest for leaves due to the high maintenance cost of rubisco (González-Zurdo et al., 2016;Suzuki et al., 2010).
For grasslands, we hypothesize that aboveground biomass mortality occurs primarily at the end of the growing season as litterfall (see next section); thus, for the canopy pools E cp = 0.95. In contrast, root mortality depends solely on turnover given that root turnover is still poorly constrained due to the complexity of factors that influence root mortality, the heterogeneity in root life spans, and difficulties in measurement techniques (Urrutia-Jalabert et al., 2017;Riley et al., 2009;Smithwick et al., 2014). Since root production, respiration, mortality, and turnover all interact, SiB4 uses root respiration efficiencies (E rp ) ranging from 0.45 to 0.6.
Studies have suggested that respiration scales predominantly with the assimilation rate (Flexas et al., 2006;Meir et al., 2008;Molchanov, 2009). SiB4 introduces an assimilation rate scaling coefficient (M A ) that is calculated separately for the canopy and root pools. For any canopy pool (cp), M A is calculated as where ClimA is the 10-year running mean of AssimD. CA Min , CA Max , CA L , and CA H are parameters. The assimilation rate scaling coefficient for the roots (M A,rp ) uses the same formulation, but with parameters specific to the root pools (RA Min , RA Max , RA L , and RA H ).
SiB4 uses two temperature processes to adjust the carbon pool turnover rates. First, freezing temperatures inhibit microbial activity slowing root respiration (Hobbie et al., 2000;Mikan et al., 2002;Schaefer et al., 2008). SiB4 exponentially decreases the freeze inhibition scaling coefficient (M F ) for any canopy pool cp using the equations M F,cp = e CF Mul (TC−CF Re ) and (21) where TC is the canopy temperature and CF Mul , CF Ref , and CF Min are parameters. This same formulation holds for the roots (M F,rp ) using root-specific parameters (RF Mul , RF Ref , and RF Min ), with M F,rp being calculated at every soil layer containing roots using the corresponding soil temperature (TD).
Second, studies have shown autotrophic respiration and root mortality increases with temperature (Bai et al., 2010;Gill and Jackson, 2000;McCormack et al., 2014). As in previous versions of SiB Denning et al., 1996;Schaefer et al., 2008), the high temperature scaling coefficient (M H ) use a Q10 relationship; however, in SiB4 M H is calculated separately for the canopy and root pools using the separate reference temperatures (CT Ref , RT Ref ), Q10 bases (CT Q10 , RT Q10 ), and maximum contributions (CT Max , RT Max ).

Litterfall
Using phenology stages updated daily provides a mechanistic approach for modeling the start of season and the progression throughout the growing season; however, changing allocation and photosynthetic rate alone combined with continuous turnover is not enough to capture the rapid grassland browning seen at the end of the growing season. To predict senescence, this approach adds three explicit processes for litterfall.
Day length, temperature, and soil moisture are commonly used environmental cues for senescence (Borchert et al., 2005;Dahlin et al., 2015;Jolly et al., 2005;Kim et al., 2015;White et al., 1997). The transfer of biomass in canopy carbon pools (cp) to dead carbon pools is calculated as where LossL is the canopy pool loss of carbon at time step t, T DL is the fraction of the canopy pool that is transferred to the dead carbon pools due to decreasing day length, T F is the fraction of the canopy pool that is transferred due to freezing temperatures, and T W is the fraction of the canopy pool that is transferred due to water deprivation.
For day length, T DL increases as the days get shorter following where LTD C , LTD R , and LTD Max are parameters. An example of daylength transfer fractions calculated at 45 • N is shown in Figure S5.
For freezing temperatures, T F increases with decreasing temperatures using a Q10 relationship following where LT Fq10 , LT Fmax , and LT Fref are parameters. An example of T F is shown in Figure S6.
For moisture deprivation, T W increases with lowering PAW such that − 1) and (28) where LTW B , LTW C , LTW R , and LTW Max are parameters. An example of T W is shown in Figure S7.

Grazing
Grasslands have a long history of grazing, and grazing is an important disturbance that has direct impacts on the carbon cycle from local-scale vegetation dynamics, to regional productivity, to global sources and sinks (Kang et al., 2007;Knapp et al., 2012;Koerner et al., 2014;Yan et al., 2013). Grazing removes aboveground plant material, which has considerable effects on the land-atmosphere interactions by changing the aboveground vegetation cover, and hence altering the LAI, litter, albedo, roughness length and resulting energy and momentum exchanges. Changing the aboveground carbon pools also alters the carbon fluxes and plant productivity, which feeds back to impact not only the root biomass, but also the resulting carbon stored in the soil. Extensive grazing occurs globally, and these local impacts scale up to have global implications.
Since grazing has widespread impacts, SiB4 includes a rudimentary grazing scheme. Following grassland management practices for sustainable grazing, the amount of removed carbon is a constant fraction of the long-term aboveground biomass (which for grasslands corresponds to the net primary production, NPP). While grazing estimates are improving, grazing intensity remains highly variable and uncertain, ranging from none to 100% in overgrazed pastures (Fetzel et al., 2017). SiB4 uses two different grazing intensities depending on climatological (10-year running mean) LAI, both of which are adjustable parameters. Currently, for grasslands with climatological maximum LAI ≥ 1, 40% of annual NPP is grazed, corresponding to moderate grazing practices. Otherwise, 10% of annual NPP is grazed. Grazing occurs daily once grasslands have a non-minimal LAI, and grazed carbon is primarily transferred to the dead carbon pools, with a small fraction immediately respired. Details for the grazing scheme are in the supplemental material (section S2).

Terrestrial Carbon Cycle Prediction
Prognostic phenology provides a powerful mechanism to integrate ecosystem processes, linking carbon pools (storage), land-atmosphere (energy) fluxes, and biosphere-atmosphere (carbon) exchanges and leading to a fully predictive terrestrial carbon cycle. Rather than relying on satellite data for the vegetation state, SiB4 determines this from fully prognostic carbon pools. Figure 2 shows an overview of the carbon pools in SiB4, along with the processes controlling their fluxes and transfers. Above and belowground live biomass and dead carbon pools can be predicted from a combination of phenophase-dependent assimilation rates and carbon allocation with environmentally dependent turnover rates, and these pools and resulting vegetation states feed back to impact land surface properties and radiative transfer, carbon and energy fluxes, and hydrology.
To combine these processes, every time step SiB4 computes the albedo, radiation, temperature, and soil moisture, along with the resulting energy exchanges, moisture fluxes, and carbon fluxes. Every day SiB4 sums the carbon fluxes and calculates carbon transfers required from disturbance (grazing). Once daily, the gains and losses are combined to update the carbon pools and phenology stage. Using the updated pools, all related land surface properties are diagnosed and used for subhourly photosynthetic assimilation as well as subhourly autotrophic and heterotrophic respiration and pool transfer rates. This sequence completes the carbon cycle, providing self-consistent predicted vegetation state, carbon pools, and land-atmosphere exchanges. An illustration of the timing sequence for the SiB4 carbon cycle along with a step-by-step outline of the calculation sequence is in section S3.

Site Demonstration
For an in-depth analysis of grassland carbon cycling using SiB4, eddy-covariance carbon fluxes, remotely-sensed vegetation states, and in situ net primary productivity data, please see Part 2 (Haynes et al., 2019). Here, we walk through one year at a single site to demonstrate the new methodology.   Knapp et al., 1998) is shown in Figure 3. During 2002, this site has two peaks in carbon fluxes and aboveground biomass to show how the phenology scheme is dynamic and capable of reacting to ecological and environmental conditions. The panels are organized in a manner to illustrate the daily calculation sequence.
Daily mean carbon fluxes (GPP, NPP, R A , and R H ) are shown in the top left panel of Figure 3. This panel shows the seasonal uptake and respiration of carbon, as well as the day-to-day variability in the fluxes. During the growth and mature leaf stages, R A is approximately half of GPP, and during the dormancy stage the maintenance costs exceed the carbon being taken up from photosynthesis due to the difficult environmental conditions.
The phenology index and the three contributing potentials are shown in the upper right panel of Figure 3. At the beginning of the year during winter, the total phenology index is near 0, since cold temperatures force the weather potential to low values. Since the LAI is minimal, the growth potential is at its maximum value; and since the days are getting longer, the day length potential is also at its maximum. In spring when the conditions become suitable for photosynthesis, the weather potential jumps up to unity, causing a maximum value in the phenology index, resulting in leaf-out. With the growth and day length indices still at their maximum values, the phenology index follows the weather potential, and the phenophase alters between leaf-out, growth, and maturity depending on daily conditions. Toward the middle of the year, the LAI becomes large enough that the growth index begins decreasing. At the same time, the site experiences heat and water stress, causing the weather potential to decrease and shifting the phenophase into dormancy. After several weeks, the environmental conditions improve, causing the weather potential to jump up to unity; however, at this

10.1029/2018MS001540
late point in the season the day length index is decreasing, which limits the phenophase to the growth stage.
As the year progresses and conditions deteriorate, both the weather potential and day length index decrease, causing the phenology to advance to the stress and dormancy stages. The last few days of the year after the solstice when the days begin getting longer, the day length index jumps back to unity to prepare for another growing season, while the phenology index remains low due to the unsuitable environmental conditions and low values of the weather potential.
The allocation fractions calculated from combining the phenology fractions with environmental adjustments are shown in the bottom right panel of Figure 3. The figure shows the pattern of leaf allocation decreasing throughout the season as more carbon is allocated to the roots. This panel illustrates the variability in carbon allocation and shows that although the main seasonal phenological progression dominates the allocation, there are contributions from the environmental adjustments, especially at the end of the growing season.
Finally, from all of the daily gains and losses, the carbon pools are updated. Sample seasonal cycles of the grassland live carbon pools are shown in the lower left panel of Figure 3. As the growing season progresses, the live carbon pools change daily. During the first half of the growing season, the above ground carbon pools grow, until the vegetation becomes mature and summer stresses begin. When this occurs, the root pool grows while the aboveground pools decrease. Once the site becomes stressed, the aboveground pools decrease both from increasing maintenance costs and from increasing carbon transfer to the dead carbon pools. In the second half of the year, the return of the growth stage from more optimal growing conditions causes a secondary peak in the aboveground carbon pools. At the end of the growing season, as the GPP continues to decrease, the costs of the aboveground live pool begin to outweigh the carbon being provided, and this combined with the high transfer rate due to stresses cause a rapid decrease in the aboveground pools.
The lower left panel of Figure  The observed and predicted phenology show similar patterns: the grassland remains dormant until mid-April, rapidly greens-up to reach a maximum in June, browns-down as the grassland experiences water and temperature stress in July, greens-up again in September due to late-season warm temperatures and rainfall in September, then finally senesces and enters dormancy in late October. However, SiB4 begins the growing season too early at this site, has too broad of a primary growing season, and senesces too rapidly at the end of the season compared to MODIS. Using the dynamic phenostage approach, SiB4 is able to simulate both the maximum LAI as well as the secondary late-season peak, in addition to the rapid transitions between these phenological events.

Conclusions
Since phenology is an important yet difficult component of models, we developed an innovative predictive phenology method that continuously predicts vegetation phenophase to link biophysical, biogeochemical, and phenological processes. The phenophase determines carbon allocation, maximum photosynthesis rate, and biomass transfer, providing a method to seasonally vary these important quantities. This dynamic approach closes the terrestrial carbon cycle using environmental mechanisms. Integrating phenology, carbon fluxes, and biomass provides a powerful constraint on the terrestrial carbon cycle, leading to an emergence of ecosystem characteristics, behaviors, and properties that can be used for evaluation, rather than relying on empirical or tuned parameters. For example, the root-shoot ratio in SiB4 is not prescribed, but instead emerges as a long-term property from daily changing carbon allocation. By viewing root-shoot ratios as emergent properties from different patterns in carbon allocation, evaluating this metric yields insight into allocation not only for different vegetation types and for different regions, but also for ecosystem responses to climate and perturbations in either time or space.
A full understanding of the linkage between phenology and terrestrial biophysics is unknown, particularly on community and ecosystem scales, thus using phenophases to link phenology with terrestrial biology provides an approach that can be used to learn about these relationships. Combining prognostic phenology with established sub-hourly photosynthesis, radiation, and hydrology processes, SiB4 is a tool that can be used to yield predictions across times scales ranging from minutes to seasons to centuries. A wide branch of research is using phenology to learn about biological mechanisms, ecological processes, and climate change, and SiB4 can be included in these studies to investigate the impact of changing environments on ecosystems, including impacts on land-atmosphere fluxes, live biomass, and carbon pools. As with any model, future work is required to improve SiB4; however, it is our hope that model improvement will coincide with learning about the processes involved in land-atmosphere interactions, the terrestrial carbon cycle, and climate change.
Evaluation at a grassland site illustrated the feasibility of this method, and a complimentary paper further evaluates this approach across grassland sites worldwide (Haynes et al., 2019). SiB4 has also been used to investigate the impact of drought in temperate grasslands across the Great Plains of North America (Curry, 2016), to quantify land surface behavior across the Sahel (Smith et al., 2018), and to investigate phenology using solar induced fluorescence (Cheeseman, 2018).