Dynamic global vegetation models underestimate net CO₂ flux mean and inter-annual variability in dryland ecosystems

We evaluated the TRENDY v7 models against net and gross CO₂ fluxes from 12 semi-arid AmeriFlux sites in the SW US that spans tree, shrub, and grass dominated sites and elevations ranging from ∼150 m to ∼3050 m (AmeriFlux Network, 2021 - figure S1A and see table S1 (available online at stacks.iop.org/ERL/16/094023/mmedia) for information on vegetation and soil characteristics for each site plus observation time period and site DOI) [24]. This part of the SW US is within the North American Monsoon region; therefore, these sites typically experience much of their rainfall during July to October, preceded by a hot, dry period in May and June ([24] figures 2(d)–(f)). The lower elevation (⩽1610 m) C3 shrub- and C4 grass-dominated sites have mean annual temperatures of 10 °C–20 °C and are predominantly driven by summer monsoon precipitation; however, winter and spring rains can also contribute to more ephemeral spring growing seasons at these sites [24, 25, 29, 30]. The high elevation (⩾1930 m) forested (conifer) sites experience cooler mean annual temperatures of <10 °C, and also have bi-modal growing seasons with available moisture coming from winter precipitation and summer monsoon rainfall [30, 31]. Sites are categorized throughout based on their mean annual NEE (see figure S1B): High elevation forest-dominated sites (US-Vcm, US-Vcp, US-Mpj, US-Fuf, US-Wjs and US-Ses) are a mean annual sink of C; whereas low elevation shrub- and grass-dominated sites (US-Wkg, US-SRG, US-Seg, US-SRM and US-Whs) 'pivot' between being a mean annual C sink or source, depending on annual water availability (figure S1B). One low elevation grassland site is a mean annual source of C to the atmosphere (US-Aud) [24].

Eddy covariance flux tower instruments at all sites collect half-hourly measurements of surface energy fluxes and NEE. NEE was partitioned into GPP and R_eco by each site PI. Gross CO₂ fluxes (GPP and R_eco) were calculated from the net CO₂ flux using the relationship between nighttime NEE and temperature [24]. Note that in this study a negative NEE implies a net CO₂ uptake into the ecosystem. Eddy covariance latent heat flux data were processed to provide evapotranspiration (ET). ET gaps were filled using a modified look-up table approach based on [32], with ET predicted from meteorological conditions within a 5 d moving window. We calculated seasonal CO₂ and water fluxes by summing the daily fluxes for the following months: November to February inclusive for the cool (winter) period; March to June hotter pre-monsoon (spring) period; and July to October for the monsoon. Note that the spring could be further split into the warm, moist months of March and April followed by the hotter, drier months of May and June; however, this entire period is marked by relatively warm and dry, moisture limited conditions compared with the winter or monsoon. To determine which seasons and gross CO₂ fluxes may be responsible for underestimate in NEE IAV, we examined the R² values obtained from the linear regression between the observed annual and seasonal C fluxes and the annual NEE.

Section 2.3.2 and tables 4 and A1 in [28, and references therein] provide details on the 14 DGVMs participating in TRENDY v7, as well as a description of the TRENDY MIP v7 protocol, including forcing data, simulation set-up, vegetation map, atmospheric CO₂ concentration data and land use change (LUC) datasets used. Simulated variables were downloaded from https://sites.exeter.acuk/trendy/. The models were forced by monthly CRU or the merged 6 hourly CRU–JRA-55 climate reanalysis datasets at 0.5 × 0.5° resolution that start in the year 1901 and have been updated to 2017 [33]. The protocol was as follows: (a) first a spinup was performed by cycling the climate forcing over the 1901–1920 period with a fixed global PFT map and atmospheric CO₂ concentration level from the year 1700 (276.59 ppm); (b) a transient simulation from 1700 to 1900 with changing atmospheric CO₂ based on proxy and measured data, transient land-use changes, and cycling of climate forcing over 1901–1920; and (c) a historical simulation with climate forcing from 1901 to 2017 with changing atmospheric CO₂, nitrogen deposition (if used) and land-use changes (S3 simulation). The spinup is run to ensure the C stocks reach equilibrium, which can be a different time period for each model but effectively allows for several thousand simulation years. Atmospheric CO₂ concentrations are based on ice-core proxy data (pre-1958) and measured atmospheric mole fraction data from the Mauna Loa and South Pole Observatory stations (post-1958) provided by the NOAA Earth System Research Laboratory [28—section 2.4.1]. The models either prescribe static PFT fractions per grid cell or simulate dynamic vegetation changes over long timescales (annual to millennial); however, historical LUC is imposed in the models. LUC is based on gross land use transitions from the land use harmonization v2 (LUH2) dataset [34] and net transitions from the HYDE (History Database of the Global Environment) v3.2 [35]. Each modeling group has its own protocol for merging this LUC data with their own PFT descriptions. Individual modeling groups also use their own expert judgment for setting their model parameter values, other required forcing data streams, and soil texture and permeable depth information. Model grid cells corresponding to each site location were extracted from the global simulations for our analysis.

Typically, model-data discrepancies are expected when comparing coarse grid (∼50 × 50 km) scale climate reanalysis forcing data used in TRENDY model simulations against in situ observations, which are representative of an area of ∼1–2 km². Therefore, we tested whether scale mismatch was responsible for model-data misfits by running the ORCHIDEE v2.0 DGVM (equivalent to 'ORCHIDEE v2.0' used in TRENDY v7) with site-based meteorological forcing, vegetation fractional cover and type, and soil texture fractions (hereafter referred to as ORCHIDEE_SL for 'site level'). Meteorological forcing data used to run the ORCHIDEE_SL simulations for each of the 12 sites included 2 m air temperature and surface pressure, precipitation, incoming long and shortwave radiation, wind speed, and specific humidity. These data were downloaded from the AmeriFlux data portal for each site (http://ameriflux.lbl.gov). The meteorological forcing data were gap-filled using downscaled and corrected ERA-Interim reanalysis data [36]. In situ measured precipitation was partitioned into rainfall and snowfall using a temperature threshold of 0 °C. We followed the same protocol for the site simulations that was used in TRENDY v7 (see section 2.3), with the following exceptions: (a) we performed an analytical 400 year spinup [37] by cycling over the available gap-filled forcing data for each site (table S1) followed by the transient and historical simulations; and (b) plant functional types (PFT) and soil texture fractions and maximum LAI were prescribed in the model based on the current data and literature for each site and did not change during the spinup, transient or historical simulations (table S1).

Several different metrics are used to evaluate the model annual NEE, GPP and R_eco, including mean bias error (MBE), the slope of the linear regression between model and observations, coefficient of determination and concordance correlation (ρ_c) [38]. Beyond our initial evaluation of the models' annual NEE using standard model evaluation metrics, we expanded our analysis in the following three additional steps: First, we examined the mean bias error (MBE) and bias contribution to the mean squared deviation (MSD—see below for MSD decomposition method) between modeled and observed NEE, GPP and R_eco to diagnose factors causing model discrepancies in mean annual NEE at the C sink and source sites. Second, we analyzed the slope of the linear relationship between model and observations, as well as the mean squared variation (MSV) contributions to the MSD (including phase and variance components—see below), to identify what is causing models' failure to capture NEE IAV. The analyses in these two steps were performed at sub-annual timescales in order to determine which seasons, and which of the gross CO₂ fluxes, are responsible for model-data discrepancies in annual NEE; thus, these two steps were intended to address objective 3 outlined in the introduction. Finally, for the key seasons controlling discrepancies NEE IAV we tested whether the models captured CO₂ flux response to variability in a range of climate drivers and measures of water availability (including observed air temperature, precipitation (rainfall and snowfall), vapor pressure deficit (VPD), soil moisture, and evapotranspiration, ET) given that variability in these drivers is a dominant control on intra- to interannual C flux variability [3, 24]. This final step was designed to address objective 4 of our study. We included an evaluation of the sensitivity of C fluxes to measured ET as a proxy of plant water availability because in water-limited ecosystems, ET is a more integrated metric of plant water availability (as opposed to soil moisture or precipitation) [24, 39] because plants have developed a range of strategies for dealing with limited water availability, including the ability to capture water from either deep or lateral root systems. Thus, ET accounts for water inputs, runoff losses and changes in soil moisture or snowpack storage, as well as differences in plant stomatal regulation and other vegetation characteristics that affect the energy balance. Another advantage of using ET is that it is measured at the same temporal scale, over the same flux footprint, and by the same instruments as NEE, whereas soil moisture is measured differently at every site.

We partitioned MSD between the model and the observations into bias, phase, and variance components [40]:

$$\begin{align}{\text{MSD}} & = \,\frac{1}{n}\mathop \sum \limits_{i = 1}^n {\left( {{x_i} - {y_i}} \right)^2} = \left( {\bar x - \bar y} \right)\nonumber\\ &\quad + \frac{1}{n}\mathop \sum \limits_{i = 1}^n {\left[ {\left( {{x_i} - \bar x} \right) - \left( {{y_i} - \bar y} \right)} \right]^2}\end{align} \tag{ 1 }$$

where x is the model estimate and y is the observations. We calculated a separate MSD for each of the TRENDY models and each variable of interest (e.g. NEE and GPP). The first term on the right-hand side of equation (1) represents the squared bias. The second term represents the MSV. The MSV represents the ability of the model to simulate variability about the mean. The MSV can be further partitioned into phase and variance components:

$$\begin{equation}{\text{MSV}} = { }{\left( {{\sigma _x} - {\sigma _y}} \right)^2} + 2{\sigma _x}{\sigma _y}\left( {1 - R} \right)\end{equation} \tag{ 2 }$$

where σ is the standard deviation and R is the correlation coefficient between the model and observations. The first term on the right-hand side in equation (2) indicates the difference between the model and observations in the magnitude of variability (i.e. the variance component, or in the magnitude of the variations) [40, 41]. The second term on the right-hand side of equation (2) represents a lack of correlation weighted by standard deviations [40]; thus, this final term can be thought of as a difference in phase between the model and observation [41]. The contributions of individual components to the overall MSD (or MSV) are calculated by dividing the bias, phase, and variance by the MSD (or MSV). We calculated all the above metrics at both annual and seasonal timescales.

Across all sites, the DGVMs generally underestimate the mean annual NEE (as seen by model values clustering around a y-axis value of 0.0 in figure 1(a) instead of on the 1:1 line, median slope values of 0 in figure 1(b), and high mean bias errors, MBE, in figure S2(a)). The models only simulate a weak C sink or source, irrespective of the strength of the observed sink/source. Similarly, the majority of models underestimate the magnitude of NEE IAV and fail to capture the correct IAV sign (as seen by slopes ≪1 in figures 1(a) and (b), green and white colors in figure S2(b), R² values ≪1 in figure S2(c) and concordance correlation, ρ_c, <0.5 in figure S2(d)).

Zoom In Zoom Out Reset image size

Decomposition of the annual NEE MSD into its component parts (see section 2.5) allows us to determine whether the MSD is dominated by the mean bias versus model inability to capture the variability about the mean (i.e. high MSV). The high elevation C sink sites have the strongest positive biases (model minus observations) in annual NEE (MBE of 100–400 g Cm⁻² yr⁻¹—figure S2(a)), while the mean C source site (US-Aud) has a negative MBE. Low elevation pivot sites have low biases because their mean annual NEE is close to zero. Biases therefore dominate the mean sink and source site contributions to MSD (bias contribution to MSD > 0.5 in figure S3(a)), whereas the pivot sites are dominated by model inability to capture variability about the mean (phase plus variance contributions to MSD > 0.5—figures S3(b) and (c)). Across all sites, models' inability to replicate the annual variability is generally more related to inaccurate phase rather than model failure to capture the magnitude of fluctuations (variance component) (cf figures S3(c) with S3(b) although there is considerable spread across models.

The site-level simulations with ORCHIDEE v2.0 (see section 2.4) yield a similar picture: ORCHIDEE_SL underestimated mean annual NEE and the magnitude of the NEE IAV at all sites (figure S4(a)). Given most of the models have similar representations of the main C cycle processes (see [28] and references therein), we suggest that these site level simulations with ORCHIDEE v2.0 demonstrate that inaccurate forcing or vegetation and soil characteristics are likely not the cause of DGVM underestimates in mean annual NEE and IAV; therefore, we need to analyze the models further to determine the root causes of these model-data discrepancies. The median slope of regression between simulated and observed annual NEE was 0.24 for both ORCHIDEE_SL and ORCHIDEE v2.0 in TRENDY v7, albeit with differences across sites, with a range from −0.11 to 1.2 for the site level simulations and −0.17–1.88 for TRENDY. Similarly, the median R² was 0.3 for both site level and ORCHIDEE_SL and ORCHIDEE TRENDY simulations, with a range from 0.01 to 0.87 for ORCHIDEE_SL and 0.004–0.8 for TRENDY. In both coarse grid and site-level simulations, the strongest biases were seen at high C sink (and source) sites (figure S4). The median slope of regression between simulated and observed annual NEE was 0.24 for both ORCHIDEE_SL and ORCHIDEE TRENDY, albeit with differences across sites (figure S4). Similarly, the median R² was 0.3 for both site level and ORCHIDEE_SL and ORCHIDEE TRENDY simulations.

To diagnose the processes that are responsible for model-data discrepancies, we examined the seasonal and gross CO₂ flux contributions to biases in the mean annual NEE. Separating annual net and gross fluxes into component seasonal fluxes (see section 2.2) allows us to identify that model underestimates of mean annual NEE at sink and source sites (positive and negative NEE MBEs, respectively) are most clearly associated with biases in spring NEE (cf figures 2(a) and (c)), with the exception of US-Vcm (strongest sink site), for which biases during the summer monsoon NEE also play a key role (figure 2(d)). Comparing the model to observed gross CO₂ fluxes we see that at high elevation sites (US-Vcm, US-Vcp, US-Mpj, US-Fuf, and, to a lesser extent, US-Wjs), the underestimate in spring NEE (positive MBEs) are mostly due to an underestimate in spring GPP (rather than an overestimate in R_eco—figure 2(c)). In the summer monsoon growing season, both GPP and R_eco are underestimated by most models (figure 2(d)), which results in monsoon NEE MBE values closer to zero (except for US-Vcm). The negative MBE in annual NEE at US-Aud (underestimate in the strength of the mean C source) across all models is associated with the negative biases in winter and spring (figures 2(b) and (c)). The models fail to capture the earlier rise in R_eco than GPP in the winter and spring at US-Aud, followed by a later increase in GPP. Instead, all models have almost the same mean seasonal trajectory in both GPP and R_eco (figure S5l). However, US-Aud was still recovering from a fire that occurred in 2002 until late 2005 [28]. Even if the models contain representations of wildfire, these schemes generally capture broad spatial fire patterns and are not likely to capture specific fires at one flux tower site; therefore, it is unlikely that the 2002 fire at US-Aud was correctly simulated, possibly explaining model discrepancies in winter and spring GPP and R_eco. Finally, we note that partitioning of NEE into the component of gross CO₂ fluxes is subject to uncertainties [42] that likely impact the absolute magnitude of the C fluxes.

Zoom In Zoom Out Reset image size

Preliminary model evaluation showed that there was no clear distinction between seasons or gross CO₂ fluxes in terms of the models' ability to capture NEE IAV: The models did not capture NEE IAV well in any season, or in either of the gross CO₂ fluxes (low slope values, and high MSV contributions to MSD across all seasons and fluxes—figure S6). Therefore, to focus our model evaluation analysis into only the seasons and gross CO₂ fluxes that are crucial for controlling NEE IAV, we performed a site-by-site analysis of the measured flux data. This analysis indicated that our model IAV evaluation efforts should focus on spring GPP at the high elevation sites, and monsoon GPP at the low elevation sites (figure S7 and see caption for more information). Models' failure to capture variability in spring GPP at high elevation forest sites (blue bars in figure 3), and variability in monsoon GPP across all low elevation grass- and shrub-dominated sites (orange bars in figure 3), are the predominant causes of their inability to replicate observed NEE IAV. High elevation sites tend to have lower slopes during the spring than the low elevation pivot sites (excluding US-Aud as a C source site—figure 3(a) cf blue bars to left of vertical grey dashed line for high elevation sites compared to those on the right for low elevation sites). Similarly, the spring GPP variance and phase components of the MSD (i.e. the MSV contribution to MSD) are generally larger at the higher elevation sink sites than at the lower elevation sites (figure 3(b)). We note that a higher MSV contribution to MSD indicates a poorer model performance in capturing variability in GPP. During the monsoon there is generally less differentiation in GPP slope values between high and low elevation sites (figure 3(a) yellow bars). Monsoon GPP slopes tend to be low (<0.5) across all site simulations, although a few exceptions exist: LPX and DLEM have GPP slopes close to 1 across at most sites, while ORCHIDEE v2.0 generally has robust GPP slope values for the low elevation sites, and LPJ for the high elevation sites (figure S8). The MSV contribution to monsoon GPP MSD is much higher (i.e. poorer performance for GPP variability) for low elevation sites (and much higher than in the spring—cf yellow bars to the right of the vertical dashed line in figure 3(b) for the low elevation sites compared to those on the left of the line for the high elevation sites). Differences in variance between the model and observations are more responsible for lack of spring GPP variability (high spring MSV) at high elevation sink sites (figure S9(a)), whereas phase differences dominate both the spring and monsoon GPP variability at low elevation pivot sites (figure S9(b)).

Zoom In Zoom Out Reset image size

To diagnose the possible cause(s) of biases in modeled spring and monsoon GPP variability, we examined the relationships between GPP and various climate drivers. Comparing the modeled and observed slope of the linear relationship between modeled GPP and observed ET (as a proxy of plant water availability—see section 2.5) reveals that all models generally underestimate the GPP response to changing plant water availability in the spring at high elevation sink sites and during the monsoon at the low elevation sites (light blue model GPP/observed ET slope ≪ red observed GPP/observed ET slope—figures 4(a) and (b), respectively). Several exceptions exist: at high elevation sites, the response of spring GPP to changing plant water availability estimated by OCN and DLEM is stronger than other models and matches the observations fairly well (figure 4(a)). The same is true for ORCHIDEE v2.0 and DLEM during the monsoon at low elevation sites (figure 4(b)). OCN had the smallest (i.e. closest to 0.0) median spring GPP MBE across all high elevation sites (followed by DLEM—results not shown). ORCHIDEE v2.0 and DLEM are two of the models that also performed well at low elevation sites in comparison to observed monsoon GPP, with model-data slopes around ∼1 (figure S8—see section 3.3). We chose to highlight the spring GPP response to changing plant water availability (ET) for the high elevation sites and monsoon GPP response for the low elevation sites in figure 4 because our analysis in section 3.3 pointed us towards these two combinations of seasons and site elevations as the most problematic in terms of simulating observed variability in GPP (also noting the importance of spring biases in GPP for capturing mean annual NEE at high elevation sites). However, in section 3.3 (figure 3) we also pointed out that the slope values between modeled and observed GPP were low for all sites during the monsoon. If we examine the GPP response to changing ET in both the spring and monsoon seasons at both high and low elevation sites (figure S10), we can see that many of the models systematically underestimate monsoon GPP moisture sensitivity at the high elevation sites (figure S10(c)), in addition to the spring. This further adds weight to what we found in section 3.3—that simulated GPP is an issue across all sites in the monsoon. Low elevation sites generally fare better in representing spring GPP moisture sensitivity (figure S10(b)), but again a few of the models are drastically underestimating this response. We also assessed the model GPP (and R_eco) relationships to other climate drivers and measures of water availability (see section 2.5). This evaluation revealed that observed sensitivity to most other climate drivers is weak or unclear (results not shown), except for the relationship between GPP and ET shown here. Therefore, we focused solely on the GPP-ET relationship. This finding was not surprising because, as discussed in section 2.5, past studies have found that ET is an integrated measure of plant water availability [24, 38]. The integrated nature of ET is further evidenced here by the fact we found a strong relationship between observed GPP and ET and only weak relationships with precipitation inputs and surface soil moisture (as well as other climate drivers such as T_air and VPD—results not shown).

Zoom In Zoom Out Reset image size