Modeling phenological responses of Inner Mongolia grassland species to regional climate change

Phenological and climate data were collected at eight (8) experimental sites in Inner Mongolia Autonomous Region of China, located between 97°10'E to 126°09'E and 37°24'N to 53°20'N (figure 1). The vegetation types at these sites include meadow steppe (found in Ergun, Ewenkizu Qi or Ewenki, and Bayaertuhushuo or Bayar), typical steppe (found in Xilinhot or Xilinhot, Xianghuang Qi or Xianghuang and Qahar Youyi Houqi or Qahar) and desert steppe (found in Wushenzhao or Wushenzhao and Alxa League or Alxa). The climate in Inner Mongolia can be characterized by cold and dry winter and warm and rainy summer, with higher rain-use efficiency than most other areas in the semi-arid and arid regions of China (Zhang et al 2011). The annual rainfall decreases gradually from the East (∼400 mm) to the West (∼200 mm) while the corresponding annual mean temperature ranges from −1.6 °C to 9.2 °C.

Zoom In Zoom Out Reset image size

The phenological data of thirteen dominant grassland species were collected at the eight experimental sites, including the timing of the main phenophases, such as green-up, heading, flowering, fruiting and senescence from 1981 to 2011 by the China Meteorological Administration and the Inner Mongolia Meteorological Bureau of China. These experimental sites were located in natural pastures with an area of ∼10 000 m² at each site where thirteen dominant grassland species (table 1) were found to include Poaceae (6 species), Rosaceae (1 species), Asteraceae (2 species), Zygophyllaceae (2 species), Fabaceae (1 species) and Tamaricaceae (1 species). A permanent experimental plot (10 m²) at each site was set up to collect phenology data for each species found within the plot. Ten individual plants for each species were examined in the afternoon every two days to record phenological information. The date of each phenological event was recorded when 50% of observed plants reached to a phenological event. For species with branches, the observations were made on the main stems only.

In addition to the phenological data, daily minimum, mean and maximum air temperatures along with daily precipitation data from 1981 through 2011 were obtained from the China Meteorological Administration recorded at nearby weather stations (supplementary table 1).

For future scenario analysis, climate data from 2041 to 2100 timeframe were taken from the output of PRECIS (Providing Regional Climates for Impacts Studies), a regional climate modeling system based on the third generation of the Hadley Centre's regional climate model (HadCM3) (Jones et al 2004). The IPCC SRES A2 scenario, which represents an environmentally pessimistic development pathway in China under medium-high emissions, is used in this study. The model predictions from 1961 to 1990 were in a good agreement with the climate record in Inner Mongolia (Xu et al 2006) and were used in estimating the grassland net primary productivity (Li et al 2014). The predicted future phenologies of the species were averaged from 2041 to 2070 to present the period of the 2050s (s indicates a decade), and from 2071 to 2100 to represent the period of the 2080s. The future phenological characteristics were compared to the historical observation data from 1981 to 2011 to derive standard deviations of phenology in the future.

In this study, the observed phenological events and climate data from 1981 to 2000 were used for model development, while the data from 2001 to 2011 were used for model verification. Since the observation data of Caragana stenophylla at Alxa was only available from 1983 to 2002, the data prior to 1995 for this species were used for model development and those acquired after 1995 were used for model validation.

In this study, we first analyzed phenological data from 1981 to 2000 using an existing thermal time (TT) model and an existing physiological time (PT) model. The PT model then was modified to include chilling effect (PTc) to assess potential impact of climate warming trend on plant phenologies. The new PT_C model was then validated using data from 2001 to 2011 and subsequently used to simulate phenologies under future climate scenarios.

The TT model is based on the concept of TT (°C d), similar to models based on growing degree days (Monteith 1977). It assumes a linear relationship between development rate and temperature above a base line T_b (Bonhomme 2000, Thompson and Clark 2006, Zhang et al 2008):

$$\begin{eqnarray}&&{\rm{TT}}=\displaystyle \sum (T-{T}_{{\rm{b}}}),\end{eqnarray} \tag{ 1 }$$

where TT is set to zero when the temperature T is below T_b.

The PT model assumes a 'broken stick' response function (see supplementary figure 1) with the temperature (Diekmann 1996, Campbell and Norman 1998, Soltani et al 2006, Zhang et al 2008, Wang et al 2014). It further assumes a maximum development rate when temperature is above a threshold and it calculates the number of days required for a certain phenological event to occur (Zhang et al 2008, Wang et al 2014). The PT can be either calculated from field observations or modeled as a function of temperature. Previous studies suggested that the function is usually a two-phase linear 'broken stick' or a smooth optimum relationship with temperature (Yin et al 1995).

It is noted that in this study we did not consider a maximum temperature tolerance (e.g. 35 °C for most plant species) that could hinder plant development, because daily temperature during the warmest month, July, in Inner Mongolia rarely reaches 23 °C. Therefore, we assume that the daily growth is at its maximum rate when the temperature is above the optimum baseline.

In this study, PT is calculated by summing daily values of the thermal effectiveness $f(T)$ (ranging from 0.0 to 1.0) starting from the first day of the year (Hanninen 1995, Cao and Moss 1997):

$$\begin{eqnarray}&&{\rm{PT}}=\displaystyle \int f(T)\cdot {\rm{d}}t,\end{eqnarray} \tag{ 2 }$$

where

$$\begin{eqnarray}&&f(T)=\{\begin{array}{c}0\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;T\leqslant {T}_{{\rm{b}}},\;\;\;\;\\ \displaystyle \frac{T-{T}_{{\rm{b}}}}{{T}_{{\rm{o}}}-{T}_{{\rm{b}}}}\;\;\;\;\;\;\;{T}_{{\rm{b}}}\lt T\leqslant {T}_{{\rm{o}}},\\ 1\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;T\gt {T}_{{\rm{o}}}.\end{array}\end{eqnarray} \tag{ 3 }$$

T is the daily temperature, T_b is the base temperature, and T_o is the optimal temperature.

In addition to the accumulation of heat needed for plants to develop leaves or flower in spring, many plants in temperate climate zone require an accumulation of cold days (chilling effect) to break plant dormancy for woody species (Campoy et al 2012, Luedeling et al 2013, Guo et al 2015) and to vernalize floral initiation in herbaceous species (Donald 1984, Cao and Moss 1997). Therefore, the PT model was further modified to account for chilling effect by assuming that lack of chilling would delay the optimal development. The modified model, PT_C, accounts for both the effect of chilling completion (CC) and thermal effectiveness to achieve phenological stages.

The chilling effect, c(T), can be modeled by setting the value to be 1.0 when the temperature is between the baseline and optimal temperatures (Cao and Moss 1997) and 0.0 when the temperature is below chilling-base temperature T_cb to account for frozen effect. When the temperature is in between T_co and T_cm, it is a simple linear function. There is no chilling effect when the temperature is above the maximum chilling temperature (T_cm) or below the chilling-base temperature (T_cb). Mathematically chilling effect, c(T), can be expressed as:

$$\begin{eqnarray}&&c(T)=\{\begin{array}{c}0.0\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;T\lt {T}_{{\rm{cb}}},\\ 1.0\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{T}_{{\rm{cb}}}\leqslant T\leqslant {T}_{{\rm{co}}},\;\;\;\;\\ \displaystyle \frac{{T}_{{\rm{cm}}}-T}{{T}_{{\rm{cm}}}-{T}_{{\rm{co}}}}\;\;\;\;\;\;\;{T}_{{\rm{co}}}\lt T\lt {T}_{{\rm{cm}}},\end{array}\end{eqnarray} \tag{ 4 }$$

where, T_cb is the base temperature to be determined experimentally, T_co is the optimal temperature, and T_cm is the maximum temperature for chilling effect.

Phenological development is dependent on CC (Cao and Moss 1997), expressed as the fraction of accumulated chilling effect c(T) from 1 October of the previous year to the total chilling requirement (CC_m). The value of CC is set to a maximum value of 1.0

$$\begin{eqnarray}&&{\rm{CC}}=\displaystyle \frac{\displaystyle \sum c(T)}{{{\rm{CC}}}_{{\rm{m}}}}.\end{eqnarray} \tag{ 5 }$$

Chilling-adjusted PT (PT_C) can be modeled as the product of daily temperature effect f(T) and CC (Cao and Moss 1997). The PT_C is the number of physiological days required under optimal development and CC. If the chilling is not completed, the development rate of species is reduced proportionally. Due to the chilling requirement of floral initiation in grass species, the chilling effect is only used for predicting flowering events. Therefore, the PT_c is expressed as:

$$\begin{eqnarray}&&{{\rm{PT}}}_{{\rm{C}}}=\displaystyle \int f(T)\times {\rm{CC}}\cdot {\rm{d}}t.\end{eqnarray} \tag{ 6 }$$

The integration time (t, days) of TT, PT and PT_C starts from 1 January of a year to each event for all species. However, the CC is calculated from 1 October of the previous year to include both initiations of leaf and the floral buds (Donald 1984, Chen et al 2014).

The base temperature (T_b) and optimal temperature (T_o) for the model were determined by fitting the model to historical observations. Doing so, we first tested T_b using a TT model (Monteith 1977, Bonhomme 2000, Thompson and Clark 2006, Zhang et al 2008), by setting T_b value from 0 °C to 10 °C to determine the T_b with the least square fitting. We then tested T_o by using the PT model at a fixed T_b determined in the first step with the value of T_o set between 10 °C and 30 °C to determine the T_o with least square fitting. To determine the threshold temperatures of chilling, we first fixed the chilling-base temperature for grass species based on previous research (Fu et al 2014), and then determined the values of optimal and maximum chilling temperatures.

Further, based on the analysis of experimental data, ${T}_{{\rm{b}}}$ and ${T}_{{\rm{o}}}$ were set to 0 °C and 15 °C for equation (3), respectively. ${T}_{{\rm{cb}}},$ ${T}_{{\rm{co}}}$ and ${T}_{{\rm{cm}}}$ were set to −5 °C, 0 °C, and 5 °C for equation (4), respectively. The same thresholds were used in the model for all herbaceous species and shrubs found at our study sites in Inner Mongolia.

The requirements of CC_m for different species were categorized into strong (CC_m = 65), moderate (CC_m = 60) and weak (CC_m = 50) sensitivities to chilling temperatures. The CC_m was estimated from the phenological data from 1981 to 2000.

The daily maximum and minimum temperatures do not change at the same rate; they are often accompanied by reduced frost occurrence and increased heat stresses (Wang et al 2015). Including the variation of minimum and maximum temperatures rather than only using daily average in the model enables better presentation of phenological characteristics in understanding climate impacts on plant growth.

The daily mean temperature was calculated as the average of the daily maximum and minimum values (Cannell and Smith 1983). The effect of diurnal temperature variation on chilling-adjusted PT (PT_C) was considered as 50% in average, 25% in both maximum and minimum air temperatures (Yang et al 2004, Zhang et al 2008). The actual daily thermal effect f(T) can be distributed among these three temperatures by:

$$\begin{eqnarray}&&f(T)=0.5\times f({T}_{{\rm{mean}}})+0.25\times f({T}_{\mathrm{max}})\\ &&\quad \quad \quad +0.25\times f({T}_{\mathrm{min}}),\end{eqnarray} \tag{ 7 }$$

where f(T_mean), f(T_max) and f(T_min) are the thermal effect calculated in equation (3), replacing T with T_mean, T_max and T_min, respectively. Calculation of actual c(T) regarding daily temperature variation is similar to the calculation of actual f(T).

To assess the predictive capability of the model, the root mean square error (RMSE) and the coefficient of variation of the root mean square error (CVRMSE) were used as indicators of model fitting:

$$\begin{eqnarray}&&{\rm{RMSE}}=\sqrt{\displaystyle \frac{\displaystyle {\sum }_{i=1}^{n}{({{\rm{obs}}}_{i}-{{\rm{pre}}}_{i})}^{2}}{n}},\end{eqnarray} \tag{ 8 }$$

where obs_i is the observed phenological date in year i at a given site, pre_i is the simulated phenological date, and n is number of years. CVRMSE (%) is calculated by the RMSE normalized to the mean of the observed values. The prediction is considered excellent if the CVRMSE < 10% and good if 10 < CVRMSE < 20% (Lin et al 2008).

In most years, the maximum temperature and precipitation occurred in July at all study sites, except at Wushenzhao station, which occurred in August. Variations in precipitation (CV = 15.0%) and mean air temperature (CV = 17.6%) across all years and sites were quite high. From 1981 to 2010, significant (P < 0.05) increasing trends were observed in air temperatures at a rate of 0.42 °C per decade for annual maximum temperature (R² = 0.30, P < 0.01), 0.47 °C per decade for the average (R² = 0.37, P < 0.01), and 0.53 °C per decade for the minimum (R² = 0.42, P < 0.01) across all study sites (figure 2(a)). In other words, minimum temperatures were increasing at a rate that is faster than mean and maximum temperatures. Annual precipitation indicate an insignificant decreasing trend (1.53 mm per decade) (R² = 0.08, P = 0.1) (figure 2(b)).

Zoom In Zoom Out Reset image size

The thermal requirement for each phenological development stage differed significantly (P < 0.01) among species (table 2). The thermal requirement for green-up ranged from 150 to 196 °C d for TT and from 9.6 to 12.4 days for PT for most species and sites, except for some species in Western regions, e.g. Cleistogenes squarrosa, Phragmites australis, Reaumuria songorica and Nitraria tangutorum, with a much higher requirement, 249–574 °C d for TT and 16.2–34.8 days for PT. The time from green-up to flowering of Stipa at five sites in middle and Eastern regions was 1575 °C d for TT and 78.8 days for PT, while that of Leymus chinensis in the same region was lower than that of Stipa, 963 °C d for TT and 52.4 days for PT. In the Western region (Alxa), the dominant species Caragana stenophylla had the earliest green-up among all species, and had the longest growing period from green-up to senescence, 3555 °C d for TT and 181 days for PT. The variation (percentage of standard error over the mean value) of physiological days for all stages was 16.9% smaller than that of TT, which indicates higher accuracy in the PT approach for quantifying grass phenology.

The maximum chilling requirements (CC_m) for 13 species in Inner Mongolia ranged from 38 to 66 (table 3), which indicates the possibility of difference in CC for different species with respect to chilling sensitivity (strong, medium or weak). The chilling time in meadow steppe or typical steppe was less than 50 while in desert steppe, it was slightly more than 60. The variation (SE/mean, %) in chilling-adjusted PT (PT_C) was 2.42% for senescence. The variations of PT_C for all developmental stages were similar to that of PT and smaller than that of TT. These results suggest that all species in Inner Mongolia meet the chilling requirement under current climatic condition.

Compared with the phenology observations from 2001 to 2011, the PT_C model results showed agreement with observations for the timing of green-up, heading, flowering, and fruiting. The overall prediction bias across all species, years, and sites was −2.4 days for time to green-up, 3.4 days for flowering, −1.7 days for fruiting, and −7.6 days for senescence. The overall CVRMSE was 9.0% for green-up, 6.9% for flowering, 6.3% for fruiting, and 6.8% for senescence (table 4). For all phenological events, the chilling-adjusted PT (PT_C) models for Stipa and Leymus chinensis had small biases (figure 3). However, the predictions at the Xianghuang site showed a large discrepancy compared to the observations (figure 3), which we postulate was caused by frequent droughts.

Zoom In Zoom Out Reset image size

Using the A2 climate scenario, we projected phenology for Stipa and Leymus chinensis in two future time periods: 2041–2070 (2050s) and 2071–2100 (2080s). For Stipa, the timing of green-up in the 2050s is 6.4 days earlier than in the 2080s, but not significantly different from historical observations. The length of the growing season of Stipa would be shortened by 9.8 days in the 2050s and by 11.0 days in the 2080s (figure 4(a)). For Leymus chinensis, the green-up time was 5.8 days earlier in the 2050s and 10.9 days earlier in the 2080s than historically observed values. The length of the growing season of Leymus chinensis is predicted to be shortened by 9.0 days in the 2050s and by 7.8 days in the 2080s (figure 4(b)). The shortening of growing season for two grass species in future climate change comparing to historical observations was due to earlier senescence.

Zoom In Zoom Out Reset image size