GEOG 566

         Advanced spatial statistics and GIScience

June 8, 2018

Interannual variation in phenology and productivity at a C3 and a C4 grassland

Filed under: 2018,Final Project @ 12:48 pm

1. Research question


Climate change is altering the production and distribution of plant species (Kelly & Goulden, 2008); however, the response of grasses to climate change is relatively understudied compared to woody plants, especially in tropical regions (e.g., Schimel et al. 2015). Globally, grasslands and savannas are estimated to comprise 30 percent of non-glacial land cover, and many important crop species are grasses (Still et al. 2003). Grasses are predicted to respond more quickly to climate change because of their short life span and dispersibility, and are aggressive invaders that can exacerbate fire cycles, cause land cover transitions, and alter carbon cycling (D’Antonio and Vitousek 1992, Angelo and Daehler 2013).

The community composition of a grassland mediates its response to its environment, and is critical to consider in forecasting the climate change impacts (Knapp et al., 2015). Grasses with the C4 photosynthetic pathway, in contrast to the ancestral C3 pathway, have comparatively higher resource-use efficienciesight, water, and nitrogen, especially under high temperatures. As a result, species using the C3 or C4 pathway will have distinct responses to a warming climate and rising CO2 (Collatz et al., 1992, 1998; Lloyd & Farquhar, 1994; Suits et al., 2005). Understanding the differential environmental constraints on the phenology, productivity, and distribution of grasses of different functional types will be critical in predicting how they will respond to future climate change, as well as identifying areas that may be prone to invasion from non-native grasses.


My initial research question, was, broadly: How are phenology and productivity of natural grasslands related to local climate variables, and how do these relationships differ between C3 and C4 grasslands?

Specifically, the research questions I answered were:

  1. Is interannual variation in production the result of multi-year patterns in production (autocorrelation)?
  2. How does the timing of the relationship between production and local climate variables differ between a C3 and a C4 grassland, and between years?

2. Data

My study sites are two eddy covariance (EC) flux tower locations in natural grassland areas located ~90 miles apart in eastern Kansas. The sites experience nearly identical climates, but the first is a natural tallgrass prairie composed of 99% C4 grass at Konza Prairie Biological Station outside Manhattan, KS, while the second is a replanted agricultural field composed of 75% C3 grass at the University of Kansas Field Station, outside Lawrence, KS. Because the two sites experience a very similar climate, despite being in distinct ecoregions, I hypothesize that photosynthetic type strongly controls differences in production and phenology at each site.

At these sites, production is measured first from satellite data using the normalized differential vegetation index (NDVI), which essentially measures the greenness of vegetation. In this analysis, I use NDVI derived from the MODIS satellites, which has 16-day temporal resolution, and is available from 2001 to present.

Production is also measured at the EC flux towers, and is calculated from ecophysiological equations that use ground-based measurements of atmospheric gas concentrations and meteorological data. The eddy covariance (EC) flux approach uses tower-mounted instruments to measure atmospheric concentrations of water and CO2, as well as air temperature, solar radiation, and other environmental data. All measurements are taken continuously every 30 minutes. EC flux data reflect the “footprint,” or area upwind of the tower where the instruments are mounted. The footprint varies with wind speed and direction, but averages about ~250m2. EC flux data span from 2008-2015.

The main metrics I am interested in are: daily total gross primary production (GPP) and the mean daily light use efficiency of production (LUE). GPP is reported in units of  gC•day^-1•m^-2. LUE is calculated as the daily sum of gross primary productivity (GPP) divided by the daily sum of the amount of incoming radiation, or photosynthetic photon flux density (PPFD). Daily LUE is converted to units of gC•MJ^-1•day^-1•m^-2 from units of µmol CO2• µmol photon^-1•day^-1• m^-2 using the molecular weight of carbon and Planck’s equation.

3. Hypotheses

My hypotheses were:

  1. Interannual variation in the phenology of production will not be strongly autocorrelated. Interannual variation in production is more strongly controlled by local climate variables
    1. a. Phenology will differ in the C3 vs C4 sites: the C4 site will have a longer growing season and a higher mean and maximum production.
  2. The C3 and the C4 site will have distinct cross-correlative relationships between production and local climate. They will differ in the timing and the strength of this relationship, and will respond distinctly to interannual variation in climate.

4. Approaches

My approach was to:

  1. Calculate annual phenology metrics from NDVI at each site, using the R package [greenbrown](
    • Phenological metrics include: the timing of the start of season and end of season, the length of the season, and mean and peak NDVI over the season. The phenological metrics are calculated by fitting curves to the annual growth cycle at each site.
    • After calculating the annual phenology metrics, I calculated autocorrelation using the acf() function n R to assess whether annual differences in phenology were a product of change over time, or cyclical trends.
  2. Calculate cross-correlation between GPP, LUE, and soil water content at each site, using the R function ccf(). Assess how cross-correlation differs between water-years  at each site.

5. Results


Fig. 1: Start, end, and length of growing season over time at each site.

Fig. 2: Mean (MGS) and Peak growing season value of NDVI at each site.

The C3 site has a consistently longer growing season than the C4 site, with an earlier start of season and an later end of season. However, based on the NDVI data, the mean and peak growing season values look like they may be the results of cyclic patterns. I fit an autocorrelation function to assess whether there were interannual patterns in the phenological indices.

Fig. 3: Autocorrelation function for start, end and length of season, at each site. Sites are plotted with distinct line types (KFS, C3 =solid; Konza, C4 = dashed).

Fig. 4: Autocorrelation function for mean and peak growing season value, at each site. Sites are plotted with distinct line types.

In the autocorrelograms above, the dashed lines represent confidence intervals for determining whether autocorrelation is statistically significant. Because the autocorrelation function rarely, and marginally, rises above the confidence intervals for any of the phenological indices, this indicates that the phenology we observe is more likely to result from interannual climatic variation– rather than cycles in production.

Production and cross-correlation with climate

After inferring that interannual variation in phenology is likely not the result of temporal autocorrelation of the indices, I investigated the relationship between production and climate. Specifically, I examined the the timing and magnitude of the relationship between GPP and soil water content. Importantly, I am comparing the timing of cross-correlation between water years, which are measured from October to September. Water years, in contrast to calendar years, are a more biologically meaningful division to use when examining the relationship between plant growth and water variables.

For this project, I will just show the relationships between GPP and soil water content. Because LUE is derived from GPP, the cross-correlation between LUE and soil water content looks similar to the relationship between SWC and GPP.

Fig. 5: Annual time series of GPP (gC / m^2/ day) and soil water content (% * 10), for 2010-2013, plotted by water year. The thick lines represent the monthly average of total daily GPP for each site, in each year. The dashed lines represent the average daily soil water content. The daily total of GPP for each site and for each year is potted in thin, light-colored lines in the background.

Inspecting the time series of GPP and soil water content, by water year, we can perceive a clear hump in SWC that precedes a hump in GPP. Calculating cross-correlation between SWC and GPP, using the ccf() function in R, allows us to quantify this relationship.

In my phenological questions, I acknowledge the caveats of NDVI data– which capture distinctions in the timing, but not the magnitude of difference between production of C3 and C4 grasses (Fig. 7). I was curious whether cross-correlation between NDVI and soil water content would reveal similar or different trends from the cross-correlation between GPP and soil water content.

Fig. 6: Time series of GPP and NDVI at the C3 and C4 sites, from 2010-2013. NDVI is multiplied by 10 in order visualize on a similar scale to GPP. Mean monthly GPP is plotted in the thick, solid lines; the daily average of GPP for each site is plotted in the background using thin, lighter-colored lines.

Importantly, NDVI does not accurately capture the magnitude of difference between C3 and C4 production. While C3 and C4 grasses have distinct relationships between NDVI and GPP, NDVI alone is not always a great indicator of production differences between the two photosynthetic types.

Fig. 7: Cross-correlation strength vs. lag for GPP and soil water content, for 2010 -2013. The cross-correlation between NDVI and soil water content for each site in each year is plotted in lighter-colored lines in the background of each plot.

The vertical lines above are plotted at x = 0, and x = 100, to simplify visual comparison. Dashed lines represent confidence intervals.

When we examine how cross-correlation between GPP and soil water content varies between sites and between years, we see a clear shift at Konza Praire, the C4 grassland in 2012. This shift indicates that GPP and LUE are peaking earlier in relation to soil water content: with a lag of about 30 days, as opposed to 100 days at the Kansas Field Station, and at Konza in previous years.

Because 2013 followed the severe drought year, 2012, this shift suggests that the C4 grasses at Konza may be more flexible in shifting the timing of their water use in response to drought or other irregular moisture conditions.

Importantly, NDVI reveals the same patterns of cross-correlation between production and soil water content.

6. Significance.

First, my results confirm that C3 and C4 grasses respond distinctly to similar climatic variation. The shift in the timing of cross-correlation of GPP with soil water content in 2013 (Fig 7) is evidence that the same climatic change will affect these photosynthetic pathways in different ways. We need to continue using both experimental and observational research to quantify the interacting effects of rising temperatures, increased aridity, and elevated CO2 on the production of C3 and C4 grasses.

Second, my results enlighten me that though NDVI does not capture differences in the magnitude of production at C3 vs C4 sites, it does capture differences in phenology / timing of production. This means that NDVI, though a coarse metric of production, is valid to use for assessing variation in the timing of production. Particularly because NDVI is widely used and easy to obtain and calculate, this distinction will be helpful in future phenology research. However, importantly, phenological metrics that calculate mean, max, or total growing season productivity– i.e., metrics that describe the magnitude of production, rather than the timing of production– should be interpreted from NDVI with caution. Ideally, other production indices should be used to assess variation across functional types in the magnitude of production.

7. My learning, w/r/t software

Over the course of this class, I have become much more confident using R Markdown to edit and generate reports for publication. This blog post was drafted and formatted using R Markdown, with only minimal modifications to fit the blog format. I have also learned:

  • how to gapfill time series data using the function zoo::na.approx, which uses linear interpolation to replace NAs in a time series
  • how to convert outputs from acf() and ccf() functions into data frames that are easily plottable using ggplot

8. My learning, w/r/t statistics

I became much more comfortable generating and interpreting correlograms. I used autocorrelation and cross-correlation to determine which variables to proceed with in my analysis.

Future directions will involve using the information I now have about cross-correlationto select variables for future regressions. E.g., I learned that GPP is cross-correlated with soil water content, and as a result, I am more interested in investigating the shifts in the timing of that relationship, than the magnitude. I learned that GPP is more strongly correlated, than cross-correlated, with air temperature and soil temperature. As a result, will use linear mixed models to assess the mean annual relationships between daily total GPP and means of these variables.

Ultimately, this work prepares me to quantify and discuss the similarites and differences in the timing of key phenological and production events for a C3 vs a C4 grassland.



Angelo, Courtney L., and Curtis C. Daehler (2013). “Upward expansion of fire‐adapted grasses along a warming tropical elevation gradient.” Ecography 36.5 , 551-559.

Collatz G, Ribas-Carbo M, Berry J (1992) Coupled Photosynthesis-Stomatal Conductance Model for Leaves of C4 Plants. Australian Journal of Plant Physiology, 19, 519.

Collatz GJ, Berry JA, Clark JS (1998) Effects of climate and atmospheric CO2 partial pressure on the global distribution of C4 grasses: Present, past, and future. Oecologia, 114, 441–454.

D’Antonio, Carla M., and Peter M. Vitousek (1992). “Biological invasions by exotic grasses, the grass/fire cycle, and global change.” Annual review of ecology and systematics 23.1 : 63-87.

Kelly AE, Goulden ML (2008) Rapid shifts in plant distribution with recent climate change. Proceedings of the National Academy of Sciences, 105, 11823–11826.

Knapp AK, Carroll CJW, Denton EM, La Pierre KJ, Collins SL, Smith MD (2015) Differential sensitivity to regional-scale drought in six central US grasslands. Oecologia, 177, 949–957.

Lloyd J, Farquhar GD (1994) 13C Discrimination during CO₂ Assimilation by the Terrestrial Biosphere. Source: Oecologia, 994, 201–215.

Schimel, David, et al. (2015) “Observing terrestrial ecosystems and the carbon cycle from space.” Global change biology 21.5 : 1762-1776.

Still CJ, Berry JA, Collatz GJ, DeFries RS (2003) Global distribution of C3 and C4 vegetation: Carbon cycle implications. Global Biogeochemical Cycles, 17, 6-1-6–14.

Suits NS, Denning AS, Berry JA, Still CJ, Kaduk J, Miller JB, Baker IT (2005) Simulation of carbon isotope discrimination of the terrestrial biosphere. Global Biogeochemical Cycles, 19, 1–15.

Print Friendly, PDF & Email


  1.   jonesju — June 15, 2018 @ 12:07 pm    

    Very good work. Next steps, possible things to consider: Need to check why the acf conf interval is +/- 0.5. Be clear to explain the difference between c3 and c4 grasses? Clarify how the sites differ with respect to C3 and c4 grasses. What is the basis for your definition of drought? How does the timing of the lag in ccf relate to the timing and direction of the soil moisture curve vs. GPP/NDVI curves? Maybe the start date of growing season is consistent between years because it is driven by radiation. Maybe take an energy budget and/or water budget perspective (including infiltration)? Look at heat energy driving the system and precipitation/ET cooling, use to explain soil moisture, and then GPP? PDSI or drought severity index? Contrast ccf of GPP with light to ccf of GPP with moisture: in those cases where curves are the same between years, can you argue that this is the driving factor because it operates consistently among years? Laura: infiltration capacity? Age of plants and rooting depth?

  2.   lundea — June 12, 2018 @ 4:31 pm    

    Hey Lila, Nice work! The information presented here contrasting NDVI with GPP seems like it will be quite valuable going forward. I was wondering, is it common to scale NDVI in order to use it as a proxy for productivity, as you did to compare curves in Figure 6? Did you have any data from 2014, so that you might look at whether the shift in lag in 2013 plays out further or only affects that year? Or do sites click back into cleaner overlapping patterns as we see in years before 2013?

RSS feed for comments on this post.

Sorry, the comment form is closed at this time.

© 2018 GEOG 566   Powered by WordPress MU    Hosted by