Leaf area index in Earth system models: how the key variable of vegetation seasonality works in climate projections

This study used the GIMMS LAI_3g v4 dataset as the observational data to evaluate the simulated LAI in ESMs. GIMMS LAI_3g is the longest available satellite LAI dataset (1982–2016), which covers the entirety of the Earth's surface with a high spatiotemporal resolution (1/12° and 15 d intervals) (Zhu et al 2013). Along with its long-term data, GIMMS LAI_3g acquires its high quality by multi-sensor-based bias corrections and 15 d maximum composites to alleviate the influence of non-vegetative effects. In addition, it substitutes possibly contaminated LAI data with snow or cloud to average the seasonal profile or for spline interpolation (Zhu et al 2013, Pinzon and Tucker 2014). Many studies have utilized GIMMS LAI_3g to understand the seasonal evolution of vegetation activity (Park et al 2018, Garonna et al 2018, Piao et al 2020a, Zhao et al 2020) and evaluate CMIP5 or land surface models (Anav et al 2013, Murray-Tortarolo et al 2013) due to its reliability and the availability of long-term data.

We used a 1 km global land cover classification dataset, generated by the Department of Geography at the University of Maryland (UMd) (Hansen et al 2000), to confine our data analysis to natural deciduous forests in the northern mid- and high-latitude regions (>30° N). Natural deciduous forests are classified into the following categories: deciduous needleleaf forest, deciduous broadleaf forest, mixed forest, and woodland. We utilized the near-surface (2 m) air temperature of the CRU TS4.01 climate dataset to identify the response of phenological dates to air temperature for the GIMMS LAI_3g dataset. CRU TS4.01 is a high-resolution (0.5° latitude–longitude grid) monthly climate dataset that covers all land areas, excluding Antarctica, from 1901 to 2016 (Harris et al 2020). It provides a credible climate dataset based on networks of surface weather station observations and is thus appropriate for examining changes in seasonal temperature.

We used the monthly LAI and near-surface air temperature outputs of a total of 14 models, with seven models participating in the CMIP5 and the remaining seven latest model versions participating in the CMIP6 (tables 1 and 2). By comparing the characteristics of the seven pairs of CMIP5 and CMIP6 models, we aimed to present the improvements in vegetation seasonality in the CMIP6 models, as compared to the previous CMIP. We analyzed near-surface air temperature and LAI from 1982 to 2014, which is the period of overlap among the CRU TS4.01 (1901–2016), GIMMS LAI_3g (1982–2016), and CMIP6 historical scenario outputs (1850–2014). Historical simulations are driven by observation-based forcing, such as greenhouse gas concentration, gridded land-use data, volcanic aerosols, sea surface temperature, and sea ice concentrations (Eyring et al 2016), thus allowing the comparison of ESM simulations with observed changes in the climate system. In the case of CMIP5, the historical scenario only spans from 1850 to 2005, which is 9 year shorter period than the CMIP6 historical scenario. To ensure the maximum overlap in the CMIP6, CRU TS4.01, and GIMMS LAI_3g outputs, we supplemented the 9 year period (2006–2014) of missing data with the output of representative concentration pathway 8.5 (RCP8.5) scenario in line with previous studies (Quetin and Swann 2018, Mudryk et al 2020). This is justified by the finding that the RCP8.5 scenario closely (within 1% for 2005–2020) tracks anthropogenic greenhouse gas emissions (Schwalm et al 2020).

Several CMIP5 and CMIP6 models provide multiple simulations based on the same experimental configuration for ensemble analyses. However, some of the models have only released their first realization, denominated as r1i1f1; thus, we only utilized the first realization, following previous research (Anav et al 2013). In order to maintain consistency in our data analysis, the GIMMS LAI_3g, UMd land cover classifications, and CMIP model outputs were regridded to a 0.5° latitude–longitude resolution grid of the CRU TS4.01 dataset. The GIMMS LAI_3g and UMd land cover were harmonized to a 0.5° resolution by averaging and majority sampling. For the LAI and temperature outputs of the CMIP models, linear interpolation was applied to transfer the relatively coarser model grid (table 1) to a finer resolution of 0.5°.

We analyzed three indicators of LAI seasonality (annual mean, amplitude, and phase) to evaluate the ESM simulation performance (figure 1). The amplitude and phase were calculated based on yearly sinusoidal components following the method of Stine et al (2009). The sinusoidal amplitude is approximately half of the LAI difference between the minimum and maximum in annual LAI cycles. The phase represents the middle of LAI seasonality, which is close to the timing of peak states, given as a degree ranging from 0° to 360°. In order to confine our study regions to forests, we excluded the regions where the annual maximum LAI was lower than 0.5 m² m⁻².

Along with seasonality indicators, annual dates for start of growing season (SOS) and end of growing season (EOS) were calculated based on the logistic curve fitting method used by Park et al (2018). This was used over the northern extratropics (>30° N) from 1982 to 2014. The SOS and EOS were computed for regions where the annual maximum LAI exceeded 0.5 m² m⁻², which is consistent with seasonality indicators. In addition, the length of the growing season (LOS) was defined as the difference between the SOS and EOS. The SOS, EOS, and LOS are phenological indicators used to understand the expansion of vegetation activity with the advancement of spring and autumn (Jeong et al 2011).

We used a Taylor diagram to show the mean spatial characteristics of seasonal indicators (annual mean, amplitude, and phase) and phenological dates (SOS, EOS, and LOS) simulated in the CMIP5 and CMIP6 models. The Taylor diagram represented both pattern correlations (as azimuthal angles) and normalized spatial standard deviations (as radial distances) to GIMMS LAI_3g. The pattern correlations to GIMMS LAI_3g were calculated from the mean patterns of each variable for each model from 1982 to 2014. The normalized standard deviations were defined as the ratio of the spatial standard deviation of each model to the spatial standard deviation of the GIMMS LAI3g. This diagram facilitates the evaluation of the simulation performance of the CMIP5 and CMIP6 models with regard to the relative spatial variability and pattern similarity. If a certain model output shows a pattern correlation and a normalized standard deviation that is close to one, the simulated distribution closely resembles the distribution of the GIMMS LAI3g with similar magnitudes of spatial variability.

We calculated T_SOS and T_EOS as the mean temperature of three months around the 33 year mean SOS and EOS, respectively, because temperature is the primary variable of vegetation seasonality in the northern extratropics (Jeong et al 2011). For the observation dataset and model outputs, we defined ΔSOS, ΔEOS, ΔT_SOS, and ΔT_EOS as the difference between the mean values within the earliest (1982–1991) and the latest (2005–2014) decade. The relationships between ΔSOS and ΔT_SOS and between ΔEOS and ΔT_EOS were then used to determine the responses of SOS and EOS to temperature variation.

To examine the mean features of simulated LAI seasonality, we compared GIMMS LAI_3g with the multi-model averaged patterns of the three seasonality indicators in the CMIP5 and CMIP6 models. Figure 2 shows the mean distributions of the annual mean, amplitude, and phase of LAI in the GIMMS LAI_3g and the arithmetic differences to those simulated by the CMIP5 and CMIP6 models (1982–2014). The spatial average of the annual mean LAI was 1.2 m² m⁻² in the GIMMS LAI_3g (figures 2(a) and (j)). However, the CMIP5 and CMIP6 models showed higher values of 2.2 and 1.7 m² m⁻², which were approximately 180% and 140% that of GIMMS LAI_3g, respectively (figures 2(d), (g) and (j)). The exaggeration of the annual mean LAI was evident in the CMIP5 models in most of the study regions, as compared to GIMMS LAI_3g. In contrast, the CMIP6 models showed a considerable improvement in terms of annual mean LAI exaggeration because the proportions of positive differences largely decreased in the study regions but remained noticeable in certain areas such as Scandinavia and Alaska. Compared to the overestimated annual mean, there was little difference between the spatially averaged amplitude in GIMMS LAI_3g (1.2 m² m⁻²) and in CMIP5 and CMIP6 (1.0 and 0.93 m² m⁻², respectively) (figure 2(k)). Regarding spatial patterns, CMIP5 and CMIP6 frequently displayed lower amplitudes along the latitude of 60° N and higher amplitudes in the subarctic regions (figures 2(e) and (h)). This indicates that the CMIP5 and CMIP6 models tend to underestimate the seasonal amplitude but exaggerate the annual mean of LAI in deciduous forests compared to GIMMS LAI_3g.

Zoom In Zoom Out Reset image size

With regard to the phase of vegetation, positive and widespread differences were shown in the study region in both CMIP5 and CMIP6 (figures 2(c), (f) and i). The spatially averaged phase in the study region was 208.8° in GIMMS LAI_3g, while the simulated phases differed by 13.1° and 13.5° in CMIP5 and CMIP6, respectively (figure 2(l)). A phase delay of 10° indicated that seasonality lags by approximately 10 d compared to the GIMMS LAI_3g. This delayed phase was clear in the boreal regions along the latitude of 60° N, and exceeded 30° in a few regions such as Alaska and Siberia, i.e. it was delayed by one month. Notably, the observed SOS has advanced by a few days over the past decades (Jeong et al 2011, Park et al 2018). However, this one-month delayed seasonality in boreal forests exhibits a significant difference in vegetation seasonality, which can distort vegetation–climate interactions during the growing season, and may mislead future projections of the impact of climate change on terrestrial ecosystems.

This delayed LAI seasonality could be a result of the inaccurate representation of the growing season, i.e. a late SOS or EOS, in the CMIP6 model. Figure 3 shows the mean patterns of SOS and EOS for 1982–2014, indicating the EOS is the primary cause of the delayed seasonality. The spatial averages of the SOSs were 152.3, 152.3, and 158.0 DOY in GIMMS LAI_3g, CMIP5, and CMIP6, respectively (figure 3(j)). In the spatial pattern displayed by the CMIP5 models (figure 3(d)), the SOS difference was relatively small, at less than 20 d in most of the high-latitude regions in Eurasia and North America. In contrast, in the CMIP6 models (figure 3(g)), the difference was larger than the results of the CMIP5, particularly in certain high-latitude regions in Eurasia. The ranges of the SOS differences from the 10th to 90th percentile were −11.4 to 11.7 DOY in the CMIP5, and −7.7 to 17.2 DOY in the CMIP6. This indicates that the CMIP6 models tend to simulate later start in the growing season dates compared to the CMIP5 with regard to the mean distribution. The majority of the CMIP5 and CMIP6 models showed relatively small biases in the mean values of SOS that were less than 20 d, regardless of the latitude (supplementary figure 3(a) (available online at stacks.iop.org/ERL/16/034027/mmedia)) and forest type (supplementary figure 4(a)). Hence, the difference in SOS is not sufficient to fully explain the delay in LAI seasonality.

Zoom In Zoom Out Reset image size

Figures 3(e) and (h) clearly show the positive differences between EOS in CMIP5 and CMIP6, which was the primary cause of the delayed LAI seasonality. The spatial averages of the EOS values were 260.8, 282.7, and 279.1 DOY in the GIMMS LAI_3g, CMIP5, and CMIP6 outputs, respectively. Positive differences in the EOS of more than 20 d were most pronounced in high-latitude regions in both Eurasia and North America, and exceeded 40 d in some regions in Alaska and Siberia. Europe and the southeastern part of North America showed negative differences in EOS, but their magnitudes and proportions were relatively minor compared to the clear positive biases shown at high latitudes. The 10th to 90th percentile ranges of the EOS difference were 5.0–34.4 d in CMIP5 and −5.0–34.0 d in CMIP6. Although the positive differences decreased slightly in CMIP6, both the CMIP models showed similar severe delays in the EOS, which led to the belated LAI seasonality simulated in the ESMs. This delay in autumn phenology, as shown in the CMIP6 models, is consistent with the results of previous studies on land surface models (Richardson et al 2012, Murray-Tortarolo et al 2013) and CMIP5 models (Anav et al 2013). Most CMIP5 and CMIP6 models showed a positive bias in mean EOS that exceeded 20 d, especially at high latitudes (>50° N; supplementary figure 3(b)). In the case of forest type, the belated EOS was distinct in deciduous needleleaf forests, mixed forests, and woodlands (supplementary figure 4(b)), which are largely distributed at latitudes higher than 50° N (supplementary figure 2). This implies that both the CMIP5 and CMIP6 models do not appropriately describe the autumn senescence of vegetation at high-latitude cold climate zones.

The delay in the EOS resulted in a relatively longer LOS in both CMIP5 and CMIP6 models compared to the GIMMS LAI_3g, with an extension of more than several pentads in certain high-latitude regions (figures 3(f) and (i)). The spatial averages of the LOS in the GIMMS LAI_3g, CMIP5, and CMIP6 were 107.8, 129.4, and 120.5 DOY, respectively. The LOS differences were mostly positive in the Northern Hemisphere in the CMIP5 outputs, and the 10th and 90th percentiles were 6.5 and 36.2 d, respectively. In the case of CMIP6, the positive differences were relatively weak compared to those of CMIP5 in the high-latitude Eurasia region, with a narrow range shown by the 10th to 90th percentile (− 4.1–28.8 d). The results reveal that the simulated phenology shown by the CMIP6 models tends to have a longer LOS resulting from a later EOS, as compared to the GIMMS LAI_3g; however, the bias to the satellite-derived LAI seasonality is decreased. The longer and later vegetation growing season contributes to the higher annual mean, weaker seasonal fluctuation, and delayed seasonality of LAI, as shown in figure 2. These large biases in the mean LAI seasonality could be caused by either a mean bias in a key climate variable (e.g. near-surface air temperature) or the misrepresentation of vegetation characteristics and processes (e.g. carbon allocation, phenology, or plant functional types). However, the CMIP5 and CMIP6 did not show a significant mean bias in near-surface air temperature, in contrast to LAI seasonality (supplementary figure 5). This suggests that the positive biases in the mean SOS and EOS were largely derived from a misrepresentation of deciduous forests in terms of vegetation properties and processes.

In addition to the spatial patterns observed in the CMIP5 and CMIP6 ensembles, we examined the LAI seasonality characteristics for each CMIP5 and CMIP6 model. Figure 4 shows a Taylor diagram based on the mean spatial patterns of seasonality indicators of phenological dates of the CMIP5 and CMIP6 models. This diagram demonstrates that the majority of the models had pattern correlations lower than 0.8, indicating that ESMs cannot accurately represent the spatial patterns of vegetation seasonality, even in historical scenarios. The number of model outputs in a relatively reliable range (with a pattern correlation of >0.6, and a normalized standard deviation of >0.5 and <1.5) totaled 20, which is approximately 25% of the total number of outputs. A comparison of the lower (30–60° N) and higher latitudes (>60° N) indicates that the CMIP5 and CMIP6 models tended to show a better simulation performance over the lower-latitude regions (supplementary figures 6(a) and (d)) in terms of the mean spatial patterns. These diagrams (figures 4(a), 6(a) and (d)) clearly illustrate that most of the ESMs faced difficulty in simulating the mean characteristics in vegetation seasonality in terms of both spatial variability and patterns.

Zoom In Zoom Out Reset image size

Figures 4(b) and (c) display colored tables containing pattern correlations of each indicator in the CMIP5 (upper panel) and CMIP6 (lower panel) models to focus on whether each model could reproduce the mean spatial patterns of LAI seasonality found in the GIMMS LAI_3g. These panels show that the ESMs simulated the SOS with the highest correlations for most of the models. Most CMIP5 and CMIP6 models showed correlations of higher than 0.5, with improved values in CMIP6. In contrast, the results of the EOS were not in agreement with the GIMMS LAI_3g in both the CMIP5 and CMIP6; most of the models showed low correlations (less than 0.7). The exception was IPSL-CM6A-LR, which showed the best simulation performance of all phenological indicators, with a credible normalized standard deviation close to one and higher pattern correlations than its previous version, IPSL-CM5A-MR. This improvement was largely due to better simulations of the EOS over the lower-latitude regions (30–60° N; supplementary figures 6(b) and (c)). However, figures 4(b) and (c) also indicates that recent models are not always better for simulating seasonality than previous versions. In the case of the multi-model ensembles, the pattern correlations were lower in CMIP6 than in CMIP5 regarding annual mean, amplitude, and EOS. This decrease in correlation was distinct for the amplitude simulated in CESM2 and NorESM2-LM, which used the community land model 5 (CLM5) as the land surface model (table 1). CESM2 and NorESM2-LM utilized the same land models and had almost identical land processes. As a result, they showed little difference in their partial correlations, despite differing evolutions of seasonal meteorology. However, compared to previous versions (CESM1-BCG and NorESM1-ME), the pattern correlations decreased significantly for most of the seasonal indicators in both CESM2 and NorESM2-LM in both higher- and lower-latitude regions (supplementary figures 6(b), (c), (e) and (f)), despite improvements in updated land surface and phenological processes (tables 1 and 2). This clearly illustrates the difficulties in vegetation modeling wherein changes in plant physiology and nutrient dynamics can result in unintended consequences in terms of vegetation seasonality.

One of the main drivers of changes in phenological dates is the near-surface air temperature for deciduous forests distributed over the northern extratropics (Jeong et al 2011). An accurate representation of responses in phenology to seasonal temperature change is closely connected to a long-term projection of changes in terrestrial ecosystems under climate change as well as its interactions with the Earth system. To briefly evaluate the phenological responses to temperature in the CMIP5 and CMIP6 models, we compared the interannual sensitivities of the SOS and EOS to the corresponding seasonal temperature (T_SOS and T_EOS) for each ESM in the deciduous forests (figure 5). The observation-based datasets (GIMMS LAI_3g and CRU TS) show a negative relationship between SOS and T_SOS, averaging −4.5 d K⁻¹ (figure 5(a)). Thus, higher spring temperatures were related to earlier SOS in deciduous forests. However, the CMIP5 and CMIP6 ensembles showed a relatively weak SOS sensitivity compared to the observation (−2.2 and −2.6 d K⁻¹ on average, respectively). In particular, BCC-CSM1-1 and NorESM1-ME exhibited the weakest sensitivities among the CMIP5 models (−1.0 and 0.0 d K⁻¹, respectively), but the newest versions indicated a stronger phenological response to changes in near-surface air temperature (−2.3 and −2.5 d K⁻¹ in BCC-CSM2-MR and NorESM2-LM, respectively). The sensitivity of SOS to T_SOS tended to be weaker in the CMIP5 and CMIP6 models, but ESMs are able to simulate the advancement of SOS due to spring warming (figure 5(a)).

Zoom In Zoom Out Reset image size

In the case of the sensitivity of EOS to T_EOS (figure 5(b)), the observed sensitivity was positive (0.5 d K⁻¹) on average, indicating that higher autumn temperatures lead to a later EOS. The EOS sensitivities to T_EOS tended to show a weaker magnitude than the SOS sensitivity to T_SOS because the EOS is affected by not only autumn temperature but also environmental conditions such as insolation, leaf aging, and water availability (Liu et al 2016). The CMIP5 and CMIP6 ensembles showed a slightly larger magnitude of the EOS sensitivity (0.7 and 0.8 d K⁻¹ on average, respectively), but the sensitivities of each model showed a large multi-model spread, from −0.3 d K⁻¹ of NorESM1-ME to 1.9 d K⁻¹ of CanESM2 in CMIP5, and from −0.6 d K⁻¹ of CESM2 to 2.9 d K⁻¹ of MIROC-ES2L in CMIP6. Notably, large proportions of deciduous forests in BCC-CSM1-1, CESM-BGC, and NorESM1-ME showed a weaker sensitivity compared to the other models, and the median sensitivities were close to zero (−0.1, 0.1, and 0.0 d K⁻¹, respectively). In particular, in the CMIP6 models, CESM2 and NorESM2-LM consistently displayed weak EOS sensitivity to autumn temperature, with medians of 0.0 d K⁻¹ for both. This irresponsive EOS in the two models seems to be related to a daylength-based trigger of growth offset in CLM4 and CLM5, and will be discussed in the subsequent sections.

Zoom In Zoom Out Reset image size

To determine whether the CMIP5 and CMIP6 models can represent recent changes in phenological dates and temperature, we examined differences in the mean distributions of SOS, EOS, T_SOS, and T_EOS between the two periods of 1982–1991 and 2005–2014 (i.e. ΔSOS, ΔEOS, ΔT_SOS, and ΔT_EOS). Figure 6(a) shows the mean and one-standard-deviation confidence ellipses of ΔSOS and ΔT_SOS (2005–2014 minus 1982–1991). The ellipse denotes the two-dimensional standard deviation of the ΔT_SOS and ΔSOS distributions; if a model displayed a change in temperature and phenological dates identical to that of CRU TS and GIMMS LAI_3g, its ellipse will be the same as the ellipse of the observation. However, only a few models overlapped considerably with the observation ellipse of ΔSOS and ΔT_SOS (figure 6(a)). MIROC-ES2L had averages of ΔSOS and ΔT_SOS closest to the observation ellipse, and a few models (CESM1-BGC, MPI-ESM-MR, and MPI-ESM1-2-LR) showed similar averages regarding the ellipse of the observation. Most models tended to exaggerate ΔT_SOS, that is, showing stronger warming around the SOS than the observation, which resulted in the distant ellipses and averages of models from the observation. In the case of ΔEOS and ΔT_EOS (figure 6(b)), five CMIP5 models (BCC-CSM2-MR, CESM1-BGC, MIROC-ESM, MPI-ESM1-2, and NorESM1-ME), and three CMIP6 models (CEMS2, MIROC-E2L, and MPI-ESM1-2-LR) exhibited means in the observation ellipse. This represents the response of the EOS to autumn warming in the models with relatively higher similarities to the observation, contrary to the inaccurate representation of mean EOS distributions (figures 3 and 4). Compared to the mean of the observation, most models tended to have a lower ΔEOS but a higher ΔT_EOS, similar to the higher ΔT_SOS shown in figure 6(a). Among the ESMs, MIROC-ES2L was shown to have the closest mean values to the observation (figures 6(a) and (b)), which indicates that MIROC-ES2L similarly simulated changes in seasonal temperature and involved phenological responses shown in CRU TS and GIMMS LAI_3g for both spring and autumn.