2.1.1 The Flexible Modeling System (FMS)

The Flexible Modeling System (FMS) is the coupler between the different model components of CM2Mc and has been developed at GFDL (Balaji, 2002).² FMS is a software framework for supporting the efficient development, construction, execution and scientific interpretation of atmospheric, oceanic and coupled climate model systems. The infrastructure is prepared to handle the data interpolation between various model grids in a parallel computing infrastructure. It standardizes the interfaces between various model components and handles the fluxes between them. The flexibility of FMS allows for the relatively simple exchange of model components. All model components are simulated on different spatial and temporal scales, and the coupler is the interface directly connected to the different parts. It interpolates the different scales to a common grid and adapts the respective fluxes to the grid of the receiving model component. Usually the variables are not directly exchanged between model components – for instance, the land model calculates the humidity of the canopy layer, and the atmosphere calculates the humidity of the lowest atmospheric layer. The coupler calculates the moisture flux between both layers and provides them to the different models on their respective spatial and temporal scales, while the different humidity variables are not transferred. By tracking these explicit fluxes of energy and water, the coupler ensures the conservation of these quantities.

2.1.2 Modular Ocean Model 5

CM2Mc employs the GFDL Modular Ocean Model (MOM) version 5 in a nominal 3^∘ × 3^∘ lateral grid, with 28 vertical levels (Galbraith et al., 2011). The meridional grid resolution increases to a maximum of 0.6^∘ at the Equator to allow the explicit simulation of some equatorial currents. The model uses rescaled pressure vertical coordinates (p^*), with the uppermost eight layers having a thickness of 10 dbar, which increases with depth to a maximum layer thickness of 506 dbar (Galbraith et al., 2011). MOM5 utilizes the tripolar model grid of Murray (1996) to avoid a singularity at the North Pole and partial bottom cells for a more accurate representation of bottom topography. Where the grid fails to resolve important exchanges of water between ocean basins, the cross-land mixing scheme of Griffies et al. (2005) is employed. MOM5 in CM2Mc is coupled to the GFDL thermodynamic–dynamic sea ice model (SIS; Delworth et al., 2006). We refer to Galbraith et al. (2011) for a more complete description of the model setup.

Enclosed in the ocean component MOM5, the Biogeochemistry with Light, Nutrients and Gases (BLING) model is run. It was developed at Princeton/GFDL as an intermediate-complexity tool to approximate marine biogeochemical cycling of key elements and their isotopes. More details can be found in Galbraith et al. (2011).

2.1.3 Atmospheric Model 2

The atmospheric module in CM2Mc is the GFDL Atmospheric Model version 2.1 (AM2; Anderson et al., 2004). It uses the finite volume dynamical core as described in Lin (2004) and implemented in CM2.1 (Delworth et al., 2006): a latitudinal resolution of 3^∘, a longitudinal resolution of 3.75^∘ and 24 vertical levels, with the lowest being at 30 m and the top at about 40 km above the surface. For the coupled setup, we use a general atmospheric time step of 1 h at which variables are exchanged with the coupler. Dynamic motion and the thermodynamic state of the atmosphere are calculated at a 9 min time step, whereas the radiation scheme has a time step of 3 h. The coupled model includes an explicit representation of the diurnal cycle of solar radiation. For a more detailed description of the model and its configuration, see Galbraith et al. (2011) and Delworth et al. (2006).

2.3.1 Interface between FMS and LPJmL5

The C main function of LPJmL5 used in the offline version is replaced by a coupler function providing the interface between the internal C functions of LPJmL5 and the Fortran functions of the CM2Mc model. The coupler function is called by FMS at an hourly time step, and it calls the specific update functions of LPJmL5 at the end of each hour, day, month or year, respectively. Ingoing and outgoing data are transferred as array arguments of this function. The mapping of the coarse resolution of the CM2Mc model to the 0.5^∘ × 0.5^∘ resolution of LPJmL5 is done by the FMS coupler. We found that the FMS land model component must be run at LPJmL5 resolution, which is 0.5^∘, so that all model components and the FMS coupler agree on which cells belong to land and which to the ocean. This yields slight changes of the land–sea mask from the original CM2M_coarse_BLING setup.

CM2Mc as well as LPJmL5 can use the message passing interface (MPI) to run the simulation in parallel on a compute cluster. CM2Mc uses FMS to set up a two-dimensional domain decomposition (i.e., it splits the global grid into rectangular domains which are mapped to concurrent MPI tasks). In contrast, the LPJmL5 grid is represented by an unsorted one-dimensional array of land cells, which is evenly distributed onto the MPI tasks. As this LPJmL5 grid is not compatible with the FMS grid exchange framework, a small wrapper library was developed for the data exchange between LPJmL5 and FMS domains. The wrapper library is called for the ingoing and outgoing data, and the time overhead for this data exchange is negligible. The coupler function as well as the wrapper library are part of the LPJmL5 distribution.

2.3.2 New canopy module

The stand-alone version of LPJmL5 does not calculate the essential coupling variables, canopy temperature and humidity, which is remedied in the coupled configuration via the addition of a new canopy module. In this new module, the canopy humidity and canopy temperature and some further quantities linked to those variables are calculated (Fig. 2). In this setting, the canopy layer corresponds to the lower boundary for the temperature in the atmosphere. The atmospheric diurnal cycle as well as the seasonal changes depend on the surface energy balance. The canopy humidity, on the other hand, is the lower boundary for the atmospheric humidity and, hence, sets the moisture content and the amount of precipitation in the atmosphere as well as the potential for evapotranspiration at the surface. A schematic overview of the different calculation steps is provided in Fig. 3.

Figure 2Schematic overview of the new canopy module.

Download

In the stand-alone version of LPJmL5, climatic input is prescribed; therefore, calculations of processes and fluxes, such as evapotranspiration, do not feed back to the atmosphere. In the coupled version, however, a small perturbation in a positive feedback loop can influence the climate and push the process towards an even larger perturbation. Thus, special attention has to be given to ensure the stability of the model by either ignoring the feedback and implementing a simple, empirical and stabilizing relationship or by increasing the complexity of the implementation, in order to get a more realistic representation of the vegetation embedded in the Earth system. The latter was done in CM2Mc-LPJmL by replacing the former simple Priestley–Taylor approach for calculating potential evapotranspiration ET₀ with the more complex and process-based Penman–Monteith evapotranspiration (Monteith, 1965). The Penman–Monteith approximation also accounts for additional parameters, such as humidity, that were previously not available in the stand-alone LPJmL5 model (Fig. 2):

where L_v is the volumetric latent heat of vaporization of 2453 MJ m⁻³, ET₀ is the evapotranspiration (in m d⁻¹), $\frac{\mathrm{d}{q}_{\mathrm{sat}}}{\mathrm{d}T}$ is the slope of the vapor pressure curve (in kPa ^∘C⁻¹), R_n is the net radiation at the surface (in MJ m⁻² d⁻¹), G is the soil heat-flux density (in MJ m⁻² d⁻¹), 86 400 is the conversion factor from seconds to daily values, ρ_a is the air density (in kg m⁻³), C_p is the specific heat of dry air ($\mathrm{1.013}\times {\mathrm{10}}^{-\mathrm{3}}$ MJ kg^−1 ∘C⁻¹), $e_{\mathrm{s}}^{\mathrm{0}}$ is the saturated water vapor pressure (in kPa), e_a is the actual water vapor pressure (in kPa), τ_aν is the bulk surface aerodynamic resistance for water vapor (in s m⁻¹) and τ_s is the canopy surface resistance (in s m⁻¹). γ is the psychrometric constant (in kPa ^∘C⁻¹) and is calculated as follows:

where P is the atmospheric pressure at the surface (in kPa), λ is the latent heat of vaporization of 2.45 MJ kg⁻¹ and μ is the ratio of molecular weight of water vapor to dry air (which is 0.622). ET₀ is presented here in the general daily form but is applied to the model on the sub-daily timescale; thus, it is divided by the number of time steps per day (in the current version 24).

Equation (1) uses the canopy surface resistance τ_s, which is the reciprocal of the non-water-stressed canopy conductance g_p (in mm s⁻¹). g_p was also slightly changed, compared with Schaphoff et al. (2018b), in order to include climate feedbacks. Following Medlyn et al. (2011), we included a PFT-specific stomatal conductance parameter g₁ (as defined in De Kauwe et al., 2015) and the vapor pressure deficit (D).

where g₀ (in mm s⁻¹) is a PFT-specific minimum canopy conductance scaled by FPC, occurring due to processes other than photosynthesis. p_a is the ambient partial pressure of CO₂ (in Pa), A_dt denotes the daily net daytime photosynthesis and 1000 is the unit conversion factor from millimeters to meters. D (in Pa) can be obtained by the canopy humidity q_ca and the saturation humidity q_sat:

While the new potential evapotranspiration is calculated at the sub-daily time step, the non-water-stressed canopy conductance is calculated at a daily time step, due to the daily calculation of the photosynthesis in LPJmL5. As climate data from FMS are available on a sub-daily basis, the photosynthesis routine uses a diurnal average of air temperature and photosynthetic active radiation.

Figure 3Schematic overview of the most important processes to determine the canopy humidity. Yellow denotes newly implemented processes in the new canopy layer in LPJmL5, green denotes internal LPJmL5 calculations and blue denotes input, provided by the FMS coupler. Daily processes are indicated by a dotted line, and processes operating at a sub-daily time step are shown using a solid line.

Download

The newly calculated potential evapotranspiration, accounting for g_p, is then also used in several LPJmL5 routines (e.g., bare soil evaporation or interception) instead of the equilibrium evapotranspiration (E_q), which was based on the Priestley–Taylor formula (Schaphoff et al., 2018b).

As a next step, we calculate the water-stressed transpiration E_tr, using the supply–demand functions of LPJmL5 as follows: the demand is calculated by the newly implemented potential evapotranspiration (Eq. 1, corrected by the fraction used for interception), and the supply is driven by vertical root distribution and phenology (as in Schaphoff et al., 2018b). The initial transpiration is then a function of the minimum of supply and demand for water. Following this, transpiration is subtracted from the various soil layers, depending on water availability. If the available water is not sufficient, transpiration decreases. The adjusted transpiration is consequently used in an inverse version of the Penman–Monteith formula in order to calculate the actual canopy conductance, linked to transpiration g_tr.

The total canopy conductance is additionally influenced by the conductance of soil evaporation (g_e) and plant interception (g_i). Therefore, we use a simple approach taking the maximum rainfall interception conductance (GI_MAX=10 mm s⁻¹) into account and considering the fraction of rainfall i stored in the canopy of a biome-dependent rainfall regime (Gerten et al., 2004):

where f_v is the vegetated grid cell fraction, and Pr is the daily precipitation. Pr and ET₀ are applied in millimeters per second (mm s⁻¹) here. The soil-evaporation conductance is calculated for the non-vegetated area of a grid cell and depends on the maximum soil conductance (GE_MAX=10 mm s⁻¹; Huntingford and Monteith, 1998) and an empirical scaling factor for the dependency of soil-evaporation conductance on soil water status (α₀=10, Zhou et al., 2006):

where w_evap is the soil water content relative to the water-holding capacity available for evaporation defined for a certain soil depth (Schaphoff et al., 2018b). Both conductances are calculated at a daily time step.

We then calculate the total canopy conductance g_c by adding g_tr, g_i and g_e and using τ_aν following Milly and Shmakin (2002).

where β_ph is the water available for photosynthesis:

with W_r as the actual soil water and $W_{\mathrm{r}}^{*}$ as the maximum available soil water. The increment of the canopy humidity q_ca per time step is then calculated as follows, using g_c:

where q_flux is the water flux from the canopy layer to the atmosphere, provided by the FMS coupler; $\frac{\mathrm{d}T}{\mathrm{d}t}$ is the gradient of the surface temperature over time; and ET is the final evapotranspiration, consisting of transpiration, evaporation, interception and sublimation from surface or vegetation into the canopy layer. For the calculation of ET, we used the Penman–Monteith equation (Eq. 1), now applying the total water-stressed canopy conductance g_c (Eq. 7). $\frac{\mathrm{d}{q}_{\mathrm{flux}}}{\mathrm{d}{q}_{\mathrm{ca}}}$ is the evaporation–humidity gradient. The total canopy conductance and the final increment of the canopy humidity, which is important for the FMS coupler, are calculated at a sub-daily time step. Equation (9) is based on Milly and Shmakin (2002) and is derived in Appendix C.

It was also necessary to implement the calculation of surface/canopy temperature within LPJmL5, which required major adaptions to the energy cycle in LPJmL5. The stand-alone LPJmL5 model calculates the temperature of different soil layers by employing a temperature transport scheme and considering air temperature as climatic input. In CM2Mc, however, the energy balance is calculated on the surface and then passed to the coupler and the atmosphere. Therefore, we had to implement this energy balance analogously in the coupled version of LPJmL5. While this surface temperature depends on several inputs from the coupler, such as radiation, it also uses several variables connected to the water cycle in LPJmL5 (evaporation, sublimation and melted water). As our approach does not account for a height-dependent canopy temperature, we used the surface temperature as an approximation for the canopy temperature, which is needed to calculate canopy humidity and evapotranspiration. Hence, surface temperature and canopy temperature are assumed to be the same, following the approach in the LaD model (Milly and Shmakin, 2002).

The soil temperature is still important for internal processes in LPJmL, such as permafrost, but it is not needed in the coupler to calculate fluxes from the land to the atmosphere. The calculation of heat transfer in the soil layers uses the heat-convection scheme as in stand-alone LPJmL5 model (Schaphoff et al., 2018b) by taking the air temperature into account, which highly depends on the canopy temperature. Both temperature calculations, for the surface/canopy temperature and for the soil temperature, operate at the fast time step.

In order to calculate the surface/canopy temperature within LPJmL5, we employed a simple energy balance formulation for the incremental change in temperature ΔT for each time step (adapted from Milly and Shmakin, 2002):

where m is the melted ice transformed to water (in kg m⁻² s⁻¹), LE_f is the latent heat of the conversion of ice into water (in J kg⁻¹), LE_v is the latent heat of the conversion of water into vapor (in J kg⁻¹), Q_sn is the energy released by snow (in W m⁻²), H is the sensible heat provided by the FMS (in W m⁻²), C_s is the heat capacity of the soil (in J kg⁻¹) and Δ_t is the fast time step duration (in s). R_n is used here in watts per square meter (W m⁻²). While the temperature is calculated individually for each stand, a weighted average over all stands within one grid cell is used in the humidity calculation and passed to the coupler. The heat balance of snow is calculated as performed for the soil layers (see Schaphoff et al., 2018b), where snow temperature changes (ΔT_snow) depend on the thermal conductivity (λ_snow=0.2 W m⁻² K⁻¹) and heat capacity (C_snow=630 000 J m⁻³ K⁻¹) of snow as follows:

The heat flux from snow (Q_snow) is calculated as follows:

where z_snow is the snow depth, T_air is the air temperature and $T_{{\mathrm{soil}}_{\left[\mathrm{0}\right]}}$ is the soil temperature of the first layer.

2.3.3 Albedo and roughness length

Albedo (β), the average reflectivity of the grid cell, is calculated as in Schaphoff et al. (2018b), based on a first implementation by Strengers et al. (2010) and later improved by considering several drivers of phenology by Forkel et al. (2014):

where the albedo for bare soil β_soil is defined as 0.3, and the albedo for snow β_snow is defined as 0.7. β_PFT is calculated for each PFT depending on the foliage projective cover (FPC) and the stem, litter and leaf albedo of the respective PFT. The value for each parameter is as in Schaphoff et al. (2018b). F_snow and F_bare are the snow coverage and the fraction of bare soil, respectively. Water bodies such as lakes and rivers have a constant albedo value of 0.1.

Roughness length z_0 m is calculated according to Strengers et al. (2010):

and

where z_b is the height of the boundary layer under stable conditions, set to 100 m (Ronda et al., 2003); $z_{\mathrm{0}\mathrm{m}}^i$ is the PFT-specific roughness length; and FPC_i is the foliage projective cover of each PFT, respectively. The coupler uses the roughness length to calculate aerodynamic resistance and surface drag and provides these variables to the different sub-models of the ESM.

2.3.4 Further changes in the coupled LPJmL5

For a global model we also need to consider Antarctica, which has not been part of the standard grid of the stand-alone LPJmL5 modeling configuration. It was implemented using a simplified approach and will be replaced with the Parallel Ice Sheet Model (PISM, Winkelmann et al., 2011) in the future. For now, Antarctica is assigned the ice soil type and a constant albedo of 0.7. The temperature balance is calculated as on the other continents.

In the stand-alone LPJmL5 model, sublimation is subsumed by a constant global value of 0.1 mm d⁻¹, likely underestimating the sublimation at high latitudes. Especially in wintertime, we do not expect much evapotranspiration; hence, the sublimation changes with meteorological conditions and becomes an important process. For this reason, we implemented the calculation of sublimation E_s using the formula from Gelfan et al. (2004):

where u is the wind speed (in m s⁻¹) from the coupler, e_s is the saturated vapor pressure (in mbar) and e_a is the air vapor pressure (in mbar).

Furthermore, the first test runs of the coupled models proved the need to tune some LPJmL5 PFT-specific parameters: we increased the effective rooting depths of the tropical tree PFTs to 2.3 m in order to counter a negative AM2 precipitation bias in northern South America. Therefore, we increased the β value of each tropical tree PFT, describing their vertical fine-root distribution in the soil column from 0.96 as in Schaphoff et al. (2018b) to 0.99 in this study.

3.1.1 Model stability

The analysis of the model stability was based on the pi-Control experiment, which ran over 750 years in total (see Sect. 2.5 for details). Here, we evaluate temperature and precipitation in terms of absolute values as well as the rate of change over time and the variability.

After the initial 300 years, the global temperature remains relatively stable at ca. 14.7 ^∘C over the remaining simulation period of 400 years with a slight drift of less than of 0.05 ^∘C per 100 years (Fig. 4a). The interannual variability in this period is ca. 0.1–0.2 ^∘C. The decreasing temperature over most of the 750-year simulation period can be explained by the energy uptake of the ocean, as deep ocean layers are not yet in equilibrium. The average precipitation follows a similar trend as temperature and reaches a relatively stable state at around 2.88 mm d⁻¹ after ca. 400 years, changing less than 0.01 mm d⁻¹ over the remaining period (Fig. 4b). The interannual variability is 0.01–0.02 mm d⁻¹.

Figure 5Yearly and decadal global mean temperature anomaly (relative to the 1951–1980 reference period) of the TR experiment of CM2Mc-LPJmL compared to GISTEMP data from 1880 to 2018. Note that, from 2004 on, only greenhouse gas forcing remains, while aerosols, solar radiation and ozone are set to their corresponding 2003 values.

Download

3.1.2 Temperature evolution over the historical period

The temperature evolution over the historical period (and hence the climate sensitivity to changes in atmospheric forcing) is evaluated by comparing the transient temperature increase in the 1880–2018 period of the TR experiment to GISTEMP evaluation data (Lenssen et al., 2019). We further evaluate the spatial impact of historical climate change without land use by comparing the years 2006–2015 of the PNV experiment with the last 10 years of the pi-Control experiment in the Supplement (Sect. S2).

The temperature evolution over the historical period from 1880 to 2018 is well captured as compared to GISTEMP evaluation data (Fig. 5). Throughout the displayed period, temperature anomalies are negative before the year 1962 and remain positive afterwards, as climate change is accelerating. While the temperature anomalies are slightly underestimated between 1980 and 2010, GISTEMP and the TR experiment both have an average global temperature increase of 0.75 ^∘C in the year 2018 relative to the reference period from 1951 to 1980. Our results are also within the range of CMIP5 models (Kattsov et al., 2013; Taylor et al., 2012, Sect. S4). The interannual variability in CM2Mc-LPJmL is ca. 0.5 ^∘C and, thus, larger than in the GISTEMP data (ca. 0.25 ^∘C), although the decadal changes are smaller in CM2Mc-LPJmL.

In the PNV experiment, climate change is also well captured, but it is weaker than when land use is included in the model (Fig. S5).

Figure 6(a) Global mean surface temperature of the TR experiment over the 1994–2003 period. (b) Surface temperature anomalies between CM2Mc-LPJmL (TR) and ERA5 data over the 1994–2003 period. (c) Latitudinal temperature mean of TR (red line) and ERA5 data (blue line) for the 1994–2003 period.

3.1.3 Surface temperature evaluation

Basic climate patterns are well captured in the annual mean surface temperature (Fig. 6a), as temperatures are increasing from polar temperatures of below −10 ^∘C towards the Equator with a maximum of ca. 25–30 ^∘C in the tropics. Desert regions are usually warmer, whereas mountainous regions are colder than the surrounding area. In the high latitudes, ocean cells are usually a bit warmer than land cells due to the ocean's ability to store heat.

Between 1994 and 2003, the average global temperature is 15.6 ^∘C compared with 14.3 ^∘C in the ERA5 data set with a NME of 0.16. While the temperatures in the tropics and temperate zone are slightly overestimated (by ca. 1 ^∘C), the poles and the boreal zone show a large negative temperature bias (up to −10 ^∘C) (Fig. 6b). The Southern Ocean has a significant positive temperature bias (ca. 3 ^∘C on average). Large differences between CM2Mc-LPJmL and ERA5 are also visible for mountainous areas, where the temperature bias is partly due to the coarse resolution of the model not adequately capturing the orographic influence of most mountain ranges on climate (e.g., the Andes or Himalaya).

While the seasonal cycle is usually well captured in CM2Mc-LPJmL, especially in Antarctica a strong seasonal temperature bias is partly balanced out in the annual mean temperature. The temperature over Antarctica is largely overestimated during the Southern Hemisphere summer, whereas it is underestimated during the Southern Hemisphere winter (Figs. S1, S2).

The latitudinal distribution of modeled mean temperature between 1994 and 2003 (Fig. 6c) shows similar values to ERA5 data from high- to midlatitudes in the Northern Hemisphere, but a slight overestimation in parts of the temperate zone and the tropics (between 70^∘ S and 40^∘ N). Specifically, the cold bias in the boreal zone leads to a slight underestimation of temperature between 60 and 90^∘ N.

The comparison of CM2Mc-LPJmL (pi-Control) and pi-CM2Mc-LaD (as in Galbraith et al., 2011) shows that similar biases in relation to ERA5 are present in both model versions. For example, both model versions slightly overestimate global temperature (Fig. S6). The strong regional biases compared with ERA5 data are also present in both model setups (Fig. S6) and are, therefore, not due to the implementation of LPJmL5.

Figure 7(a) Global mean precipitation of the TR experiment for the 1994–2003 period. (b) Precipitation anomalies between CM2Mc-LPJmL (TR) and ERA5 data over the 1994–2003 period. (c) Latitudinal temperature mean of TR (red line) and ERA5 data (blue line) over the 1994–2003 period.

3.1.4 Precipitation evaluation

The spatiotemporal pattern of global precipitation is well simulated with a global average of 2.86 mm d⁻¹ and a maximum of up to 10 mm d⁻¹ in the tropics close to the Intertropical Convergence Zone (ITCZ; Fig. 7a). Regions with little to no vegetation, such as deserts and polar areas, receive very little precipitation throughout the year.

Precipitation biases with respect to ERA5 data are, however, stronger than temperature biases, with an NME of 0.50 compared with 0.16 for temperature (Fig. 7b). The biases are strongest at the Equator with an apparent shift of the ITCZ. While precipitation in the Pacific is underestimated directly at the Equator, it is overestimated north and south of the Equator (Fig. 7b). Moreover, northern South America shows a large negative precipitation bias.

The seasonal patterns (Figs. S3, S4) confirm the imprecise modeling of the ITCZ, which remains north and south of the Equator for a large part of the year and passes the Equator region relatively swiftly. While precipitation south of the Equator is overestimated, it is underestimated north of it.

The latitudinal annual mean precipitation between 1994 and 2003 (Fig. 7c) compares well with observations, displaying the global precipitation maximum in the tropics, local minima in the subtropics and very low values at high latitudes. The tropics, however, show a shifted maximum. While the ERA5 global precipitation maximum over the Pacific is at ca. 10^∘ N and a local smaller maximum is at −10^∘ S, CM2Mc-LPJmL models the global maximum at roughly −10^∘ S and a smaller local maximum at ca. 10^∘ N. The difference between the two maxima is less pronounced compared with ERA5.

The comparison of the results of CM2Mc-LPJmL with the original pi-CM2Mc-LaD model shows similar biases in relation to ERA5 for both model versions. Neither of the models precisely captures the behavior of the ITCZ, especially over the Pacific. Both models also show a large dry bias in northern South America (Fig. S6).

Figure 8(a) Mean global aboveground biomass of GlobBiomass evaluation data. (b) Mean global aboveground biomass of CM2Mc-LPJmL (TR) over the 2006–2015 period. (c) Difference in the aboveground biomass between CM2Mc-LPJmL and GlobBiomass evaluation data. Blue (red) denotes an overestimation (underestimation) of biomass by CM2Mc-LPJmL. (d) Latitudinal sum of aboveground biomass from CM2Mc-LPJmL (blue line, R²=0.64, NME = 0.56), stand-alone LPJmL5 (black line, R²=0.94, NME = 0.35) input data and GlobBiomass evaluation data (red line).

3.1.5 Vegetation cover and biomass

While the evaluation of temperature and precipitation is performed for the years 1994–2003, we compare average model results for aboveground biomass (AGB) and the dominant PFT for the years 2006–2015 due to availability of evaluation data.

Simulated AGB shows a good pattern overall, with the largest values in the tropics, decreasing biomass in the subtropics, and a local maximum in the temperate and boreal zone (Fig. 8d). In vegetation-free areas, such as deserts or polar regions, simulated AGB is zero or very close to zero (less than 200 g C m⁻²). When comparing AGB against GlobBiomass (Fig. 8a), spatial differences emerge (Fig. 8c). While simulated AGB is slightly overestimated in boreal North America and Asia, it is underestimated in the European temperate zone and in Scandinavia, extending into eastern Europe and West Siberia. In most of the other temperate, Mediterranean-type and subtropical regions, AGB matches the observed values. In the tropics, AGB is overestimated in semiarid regions, whereas wet tropical rainforests are mostly underestimated, especially the eastern Amazon. AGB shows good agreement in the seasonally dry Cerrado region in South America but appears to be overestimated in the Caatinga in northeastern Brazil. In central Australia, AGB matches observations but is overestimated in the northern part of the continent and underestimated in the southeastern part (Fig. 8c).

Figure 8d compares the latitudinal mean of CM2Mc-LPJmL and LPJmL-offline with the evaluation data. LPJmL-offline shows better performance than the coupled model, with a smaller NME (0.35 vs. 0.56) and a better R² (0.94 vs. 0.64). While both models underestimate biomass in the tropics, biomass in the boreal zone is overestimated by CM2Mc-LPJmL and underestimated by the stand-alone LPJmL5 model compared with GlobBiomass. The LPJmL5 stand-alone version is forced by a reanalysis climatic input with a 0.5^∘ spatial resolution, and the model is calibrated to this specific climate conditions; therefore, better model performance is expected. Modeled biomass is also within the range of CMIP5 models (Fig. S7).

Figure 9(a) The dominant PFT for each cell, modeled by CM2Mc-LPJmL. Cells with more than 50 % land use are masked using gray. Cells with less than 200 g C m⁻² are shown using white. Full names of PFTs can be found in Appendix B. (b) The sum of squares errors of ESA CCI land cover for each PFT in each cell. Blue areas have a small error, and red areas have a large error. The error shown here is absolute; hence, areas with a low PFT cover for both model and evaluation data are small compared with areas with a large PFT cover.

The geographic distribution of dominant PFT cover in CM2Mc-LPJmL follows the spatial pattern of the biomass distribution (Fig. 9a). The tropics are mostly dominated by the evergreen tree PFT. In the tropical savanna areas, the tropical deciduous tree PFT dominates, along with the C₄-grass PFT. The temperate zone is dominated by land use with some summergreen trees – most common in areas such as Europe. The boreal zone is correctly covered by boreal needle-leaved and boreal summergreen trees, and the tundra zone is covered with polar grasses. To better visualize the model error for the PFT distribution, we produced an error map that consists of the sum of the square error for each PFT per cell (Fig. 9b). In tropical rainforests, the error with respect to the evaluation data is relatively small. Drier savanna areas show a much larger error, as do parts of the temperate and the boreal zone. Areas with a small FPC fraction show a small error because the error metric takes absolute errors into account. This applies to desert regions in Africa, the Arabian Peninsula and central Australia.

Figure 10Difference between the LU-only (2006–2015) and the pi-Control (last 10 years) experiments for (a) mean precipitation, (b) mean surface air temperature and (c) mean aboveground biomass.