LPJmL4 – a dynamic global vegetation model with managed land – Part 2: Model evaluation

Abstract

The dynamic global vegetation model LPJmL4 is a process-based model that simulates climate and land use change impacts on the terrestrial biosphere, agricultural production, and the water and carbon cycle. Different versions of the model have been developed and applied to evaluate the role of natural and managed ecosystems in the Earth system and the potential impacts of global environmental change. A comprehensive model description of the new model version, LPJmL4, is provided in a companion paper. Here, we provide a full picture of the model performance, going beyond standard benchmark procedures and give hints on the strengths and shortcomings of the model to identify the need for further model improvement. Specifically, we evaluate LPJmL4 against various datasets from in situ measurement sites, satellite observations, and agricultural yield statistics. We apply a range of metrics to evaluate the quality of the model to simulate stocks and flows of carbon and water in natural and managed ecosystems at different temporal and spatial scales. We show that an advanced phenology scheme improves the simulation of seasonal fluctuations in the atmospheric CO${}_{\mathrm{2}}$ concentration, while the permafrost scheme improves estimates of carbon stocks. The full LPJmL4 code including the new developments will be supplied open source throughhttps://gitlab.pik-potsdam.de/lpjml/LPJmL. We hope that this will lead to new model developments and applications that improve the model performance and possibly build up a new understanding of the terrestrial biosphere.

Introduction

The terrestrial biosphere is a central element in the Earth system supporting ecosystem functioning and also providing food to human societies. Dynamic global vegetation models (DGVMs) have been developed and used to study biosphere dynamics under climate and land use change. LPJmL4 is a DGVM with managed land that has been developed to investigate the potential impacts of climate change on the terrestrial biosphere, including natural and managed ecosystems, and is now described in full detail in the companion paper . LPJmL and its predecessors were originally benchmarked against ecosystem carbon and water fluxes and global maps of vegetation distribution , run-off , agricultural yield statistics , satellite observations of fire activity , permafrost distribution and active layer thickness , satellite observations of fraction of absorbed photosynthetically active radiation (FAPAR) and albedo , and atmospheric CO${}_{\mathrm{2}}$ concentrations . These previous evaluation studies focussed on single processes or components of the model. Here we present a comprehensive multi-sectoral evaluation to demonstrate that LPJmL4 can consistently represent multiple aspects of biosphere dynamics.

LPJmL4 spans a wide range of processes (from biogeochemical to ecological aspects, from leaf-level photosynthesis to biome composition) and combines natural ecosystems, terrestrial water cycling, and managed ecosystems in one consistent framework. As such, it is increasingly applied for cross-sectoral studies, such as the quantification of planetary boundaries and SDG interactions , and the multidimensional impacts of climate and land use change e.g.. With this complexity, its evaluation against historical observations along multiple dimensions is essential . For such a purpose, standardized benchmarking systems have been proposed and iLAMB (https://www.ilamb.org/), the international land model benchmarking project, has been established. In the present evaluation of a broad range of fundamental features of the LPJmL4 model, we basically follow the benchmarking procedures, variables, performance metrics, and diagnostic plots suggested by and . Thus the presented evaluation goes well beyond earlier evaluations of DGVMs and LPJmL (and its predecessors) itself. We pay special attention to LPJmL4's capability to reproduce observed seasonal and inter-annual dynamics and patterns of key biogeochemical, hydrological, and agricultural processes at various spatial scales. In doing so, we highlight the model's unique feature of representing the interaction of processes for both natural and agricultural ecosystems in a single, internally consistent framework.

Model benchmark

In the following we describe in detail the model benchmarking scheme employed here, which allows for a consistent evaluation of processes simulated by LPJmL4 at seasonal and annual resolution and at spatial scales from site level (using e.g. eddy flux measurements for comparison) to global level (using e.g. remote sensing products). The evaluation spans the time period from 1901 to 2011. The benchmarking analysis also considers results from different model set-ups and previous model versions in order to demonstrate advancements achieved with the current LPJmL4 version and the sensitivity of results to individual new modules.

Model set-up and simulation experiments

As described in , we drive the model simulations with observation-based monthly input data on daily mean temperatures from the Climatic Research Unit (CRU TS version 3.23, ) and precipitation provided by the Global Precipitation Climatology Centre (GPCC Full Data Reanalysis version 7.0; ). Shortwave downward radiation and net downward longwave radiation are reanalysis data from ERA-Interim . Monthly average wind speeds are based on the National Centers for Environmental Prediction (NCEP) reanalysis data and were regridded to CRU (NOAA-CIRES Climate Diagnostics Center, Boulder, Colorado, USA; ). The number of wet days per month, which is used to allocate monthly precipitation data to individual days of the corresponding months, is derived synthetically as suggested by . Dew-point temperature is approximated from daily minimum temperature . Global annual values for atmospheric carbon dioxide concentration are taken from the Mauna Loa station . The spatial resolution of all input data is 0.5${}^{\circ }$ and the model simulations are conducted at this spatial resolution. All model simulations are based on a 5000-year spin-up simulation after initializing all pools to zero. A second spin-up simulation of 390 years is conducted, in which human land use is introduced in 1700, using the data of . In addition to the original dataset description of , sugar cane is now represented explicitly. Cropping intensity as calibrated following is kept static in the simulations, whereas sowing dates are computed dynamically as a function of climatic conditions until 1971 following and kept static afterwards. Soil texture is given by the Harmonized World Soil Database (HWSD) version 1 and parameterized based on the relationships between texture and hydraulic properties from . The river-routing scheme is from the simulated Topological Network (STN-30) drainage direction map . Reservoir parameters are taken from , and locations are obtained from the GRanD database . We test the influence of specific processes that have been implemented or improved in LPJmL4 (specifically, permafrost, phenology, and fire) on overall model performance by conducting the following factorial experiments.

• LPJmL4-GSI-GlobFIRM is a simulation with all standard model features enabled as used in , i.e. with land use, permafrost dynamics, the growing season index (GSI) phenology scheme, and the simplified fire model (GlobFIRM). This model experiment is the default LPJmL4 model experiment.

• LPJmL4-GSI-GlobFIRE-PNV is the same, but for potential natural vegetation (PNV) to evaluate the role of managed land in global pattern and processes. This model experiment mimics the original LPJ model (i.e. without agriculture) but with improved phenology.

• LPJmL4-NOGSI-GlobFIRM is a simulation with land use, permafrost dynamics, and the simplified fire model, but without the GSI phenology for testing the sole effect of the GSI phenology. Instead of the GSI phenology, here we use the original phenology model that is based on a growing-degree day approach. This experiment mimics the LPJmL 3.5 version (including the LPJ core, agriculture, and permafrost) as described in .

• LPJmL4-NOGSI-NOPERM-GlobFIRM is a simulation with land use and the simplified fire model but without permafrost and without the GSI phenology. This model experiment mimics the original LPJmL 3.0 model with the LPJ core and the agricultural modules .

• LPJmL4-GSI-SPITFIRE is a simulation set-up as LPJmL4-GSI-GlobFIRM but with the process-based fire model (SPITFIRE; ). This experiment is an LPJmL4 model run with an alternative fire module.

Evaluation datasets

Following we compare LPJmL4 simulations against independent data for vegetation cover, atmospheric CO${}_{\mathrm{2}}$ concentrations, carbon stocks and fluxes, fractional burnt area, river discharge, and FAPAR. Beyond these suggestions of , we extend the benchmarking system to datasets of eddy flux tower measurements of evapotranspiration and net ecosystem exchange rate (NEE). Ecosystem respiration ($R_{\text{e}}$) is evaluated against both eddy flux measurements and operational remote sensing data. Crop yields are evaluated against FAOSTAT data . For FAPAR, we use not just one but three different reference datasets to account for uncertainties from multiple satellite datasets (see Sect. ). We also compare LPJmL4 results against data that are not fully independent of other models (mostly empirical, data-driven modelling concepts), acknowledging the limitations of these data in a benchmark system. However, this allows for the assessment of LPJmL4's performance in additional aspects for which fully data-based products are not available. These data comprise global gridded datasets of vegetation or aboveground biomass carbon , cropping calendars , global gross primary production (GPP) , $R_{\text{e}}$ , soil carbon , and evapotranspiration . We use both site-level and global gridded data because they provide complementary information but have different advantages for the comparison with simulated data like those from LPJmL4. Site-level data are fully independent from model estimates and assumptions, but typically only represent a specific ecosystem with a certain vegetation and soil type and a specific site history. Thus site-level data have only a limited representativeness for 0.5${}^{\circ }$ grid cells. On the other hand, global gridded data of GPP and $R_{\text{e}}$ are available at the same scale and thus can be directly compared to simulation outputs of DGVMs. However, global gridded datasets usually rely on empirical modelling approaches and ancillary data to upscale and extrapolate site-level data to large regions. Nevertheless, specific site conditions like forest management affecting site age, biomass, and carbon fluxes can hardly be re-simulated for a large number of global sites within a DGVM. Although reject the use of such datasets for model benchmarking because they depend on modelling approaches, we accept the additional use of such datasets because they prevent the scale mismatch between site-level data and global DGVM simulations.

Vegetation cover

We compare simulated vegetation cover to the ISLSCP II vegetation continuous fields of as suggested by . This dataset is a gridded snapshot of vegetation cover for the years 1992–1993 from remote sensing data and distinguishes bare soil, herbaceous, and tree cover fractions aggregated to 0.5${}^{\circ }$ resolution . Tree cover fractions are further distinguished into evergreen vs. deciduous and into broadleaved vs. needle-leaved tree types, respectively. The herbaceous vegetation class includes woody vegetation that is less than 5 m tall. Data uncertainties increase in regions where tree cover is $<$ 20 $\mathrm{\%}$ due to understorey vegetation and soil disturbing the signal, as well as above 80 $\mathrm{\%}$ due to signal saturation . To test if the simulated land cover of LPJmL4 performs better than a randomly generated land cover distribution we compare the performance of LPJmL4 to the random model as suggested by Sect. 2.3.5, whereas in the original dataset ISLSCP II vegetation continuous fields were randomly resampled.

Atmospheric CO${}_{\mathrm{2}}$ concentration

To evaluate the model's capacity to capture global-scale, intra- and inter-annual fluctuations of atmospheric CO${}_{\mathrm{2}}$ concentrations as driven by the uptake activity of the terrestrial biosphere, we compare simulated CO${}_{\mathrm{2}}$ concentrations with those recorded continuously at two remote measurements at Mauna Loa (MLO; 19.53${}^{\circ }$ N, 155.58${}^{\circ }$ W) and Point Barrow (BRW; 71.32${}^{\circ }$ N, 156.60${}^{\circ }$ W; see for further details on these measurements). We use monthly CO${}_{\mathrm{2}}$ concentrations from flasks and continuous measurements from 1980 to 2010 for the comparison with LPJmL4 simulations. CO${}_{\mathrm{2}}$ observations were temporally smoothed and interpolated using a standard method . The atmospheric transport model (TM3; ) in Jacobian representation simulates the global CO${}_{\mathrm{2}}$ transport using estimates of net biome production (NBP; here simulated by LPJmL4; see ), estimated net ocean CO${}_{\mathrm{2}}$ fluxes from the Global Carbon Project and fossil fuel emissions from the Carbon Dioxide Information Analysis Center (CDIAC; ). Atmospheric transport in TM3 is driven by wind fields of the NCEP reanalysis at a spatial resolution of $\mathrm{4}{}^{\circ }\times \mathrm{5}{}^{\circ }$.

Terrestrial carbon stocks and fluxes

Model-independent reference data for carbon stocks and fluxes are available from for various sites globally distributed. This dataset is comprised of vegetation carbon, aboveground biomass, GPP, and net primary production (NPP). GPP flux data from are based on eddy flux measurements and are subject to the uncertainties reported in Table 2. Contrastingly, NPP data are derived from direct measurements of continuous leaf-litter collection, allometry-based estimates of stem and branch NPP from basal measurements, root NPP estimates from soil cores, mini-rhizotrons or soil respiration, and destructive understorey harvest. Estimates here are subject to uncertainties, depending on the sampling methods . Several individual sites of this dataset can be located within one simulation unit of a 0.5${}^{\circ }$ grid cell and we thus compare simulated values to the range of site measurements in that grid cell. Alternatively to the site-based GPP data from , we also compare spatial patterns and grid-cell-specific GPP simulations to the GPP dataset of , as also suggested by . This global dataset is based on a larger set of eddy flux tower measurements than the dataset of , but uses additional satellite and climate data and empirical modelling for extrapolation to full global coverage. $R_{\text{e}}$ is evaluated for the time period 2000 to 2009 directly against plot-scale FLUXNET (http://fluxnet.fluxdata.org/data/la-thuile-dataset/) measurements , but also against large-scale $R_{\text{e}}$ estimates from an empirical model based on operational remote sensing data by the Moderate Resolution Imaging Spectroradiometer (MODIS) with a resolution of 1 km and 8 days . In addition to GPP, $R_{\text{e}}$, and NPP, we also compare simulated NEE fluxes with eddy flux tower measurements directly. We use 70 time series of estimated NEE from eddy flux tower sites that measure the exchanges of carbon and water fluxes continuously over a broad range of climate and biome types . Nevertheless, eddy flux tower sites are not well distributed across the globe and sites in the temperate and boreal zone are better represented than the tropical zone. For the global comparison of the soil and vegetation carbon stocks we use the data compiled by . The soil organic carbon (SOC) estimations are based on the Harmonized World Soil Database (HWSD) . used an empirical model to calculate SOC stocks ($\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$) from soil organic content ($\mathrm{\%}$), layer thickness ($\mathrm{m}$; here for the first 3 $\mathrm{m}$), gravel content (% vol), and bulk density ($\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$). They pointed out that regions such as North America and northern Eurasia are less reliable as HWSD was a work in progress at that time. The vegetation carbon data of are based on a forest biomass map for temperate and boreal forests from microwave satellite observations , a biomass map for tropical forests based on lidar observations , and an additional estimate of grassland biomass. Uncertainties in biomass are in most regions between 30 and 40 $\mathrm{\%}$ and are strongly related to uncertainties in belowground biomass. We also compare simulated aboveground biomass to the estimates of , which is also based on satellite-based passive microwave data. This comparison requires additional assumptions on the separation of aboveground and belowground biomass in LPJmL4 simulations. estimate for 2000 a global aboveground biomass of 362 $\mathrm{Pg}\mathrm{C}$ with a 90 $\mathrm{\%}$ confidence interval of 310–422 $\mathrm{Pg}\mathrm{C}$.

Terrestrial water fluxes

River discharge measurements are taken from theArcticNET (http://www.r-arcticnet.sr.unh.edu/v4.0/index.html) and UNH/GRDC (http://www.grdc.sr.unh.edu/index.html) datasets for 287 gauges . From this database, we only selected river gauges with catchment areas $\ge$ 10 000 ${\mathrm{km}}^{\mathrm{2}}$ as the model set-up and resolution are not suitable for comparison with smaller catchments. We also only selected river gauge records with a temporal coverage of more than 95 $\mathrm{\%}$ of the observation period and an observation period longer than 2 years at a monthly resolution.

Evapotranspiration fluxes are taken from the FLUXNET database (http://fluxnet.fluxdata.org/data/la-thuile-dataset/) and comprise 126 sites, of which we selected sites ($n=\mathrm{99}$) with at least 3 years of data available . Additional to site-level data, we used global gridded ET data from , which is based on an upscaling of site-level eddy covariance observations with satellite and climate data using a machine learning approach.

The irrigation withdrawal and consumption data that we compare to are from other modelling approaches. Nonetheless, human water use for irrigation is an important component in the terrestrial water cycle and we discuss modelled LPJmL4 estimates in comparison to other model-based estimates, acknowledging the limitation of this comparison and addressing different sources of uncertainty.

Permafrost

For the evaluation of simulated permafrost dynamics, we use the measured thaw depth data from 131 stations of the Circumpolar Active Layer Monitoring (CALM) station dataset (https://www2.gwu.edu/~calm/ ) and the International Permafrost Association (IPA) Circum-Arctic Map of Permafrost (http://nsidc.org/data/ggd318; ). The distribution of permafrost is based on regional elevation, physiography, and surface geology. The permafrost extent represents four classes which categorize the percentage of the ground underlain by permafrost (continuous, 90–100 $\mathrm{\%}$; discontinuous, 50–90 $\mathrm{\%}$; sporadic, 10–50 $\mathrm{\%}$; isolated patches of permafrost, 0–10 $\mathrm{\%}$).

Fractional area burnt

For the evaluation of simulated fire dynamics, we employ data on fractional area burnt from the Global Fire Emissions Database GFED4 version 4 (GFED4; http://www.globalfiredata.org/; ) for the period 1995 to 2014 and Climate Change Initiative (CCI) Fire version 4.1 (http://cci.esa.int/data; ) for the period 2005 to 2011. Mean annual burnt area was computed for both datasets for the overlapping period (2005–2011). Both datasets are derived from satellite data. Active fire data were used in GFED4 to prolong the dataset prior to the MODIS period (i.e. for 1995–2000).

Fraction of absorbed photosynthetic active radiation and albedo

Data on the fraction of absorbed photosynthetically active radiation (FAPAR) are derived from three different satellite datasets to account for differences between datasets for model evaluation (see Table , ): the MODIS (USGS, 2001) FAPAR , the Geoland2 BioPar (GEOV1) FAPAR dataset (hereafter called VGT2 FAPAR), and the GIMMS3g FAPAR dataset . The MODIS FAPAR dataset is taken from the MOD15A2 product with a temporal resolution of 8 days at a spatial resolution of 1 km, covering the period 2001 to 2011. VGT2 is based on SPOT VGT with a temporal resolution of 10 days and 0.05${}^{\circ }$ spatial resolution , covering the period 2003 to 2011. The GIMMS3g dataset has a 15-day temporal resolution and $\mathrm{1}/\mathrm{12}$${}^{\circ }$ spatial resolution and covers the period from 1982 to 2011. Data on FAPAR are also subject to uncertainties from the processing of the remotely sensed data and are not available continuously for all areas. We compare the spatial patterns of the peak FAPAR, the temporal dynamics of FAPAR in each grid cell, and seasonal variations in FAPAR averaged for Köppen–Geiger climate zones for the three different FAPAR datasets. The aggregated FAPAR represents the average monthly time series for all grid cells that belong to a certain Köppen–Geiger climate zone (see also ). For the Köppen–Geiger climate zones, FAPAR time series are averaged over all grid cells that belong to the same Köppen–Geiger climate zone (see also ). For the evaluation of the reflectance of the Earth's surface we used the MODIS C5 albedo time series dataset (https://lpdaac.usgs.gov/dataset_discovery/modis/modis_products_table/mcd43c3) from 2000–2010 , which we also aggregated to Köppen–Geiger climate zones for the evaluation here.

Agricultural productivity

Detailed data on crop growth and productivity are available for individual sentinel sites . For global-scale or regional simulations, reference data are available only for crop yields and in (sub-)national aggregations e.g. or as processed and interpolated gridded products . In all yield data statistics outside of well-controlled field experiments, yield levels and inter-annual variability are not only affected by variability in weather, but also by variance in management conditions, such as sowing dates, variety choices, cropping areas, fertilizer inputs, and pest control . Consequently, it is difficult to evaluate model performance from a comparison of simulated yields with static assumptions on most management aspects with yield statistics in which the contribution of weather variability to yield variability is unknown. propose a combination of global gridded crop model simulations and different observation-based yield datasets to establish a benchmark for global crop model evaluation. Generally, global gridded crop models perform well in most regions for which statistical models can detect a significant influence of weather on crop yield variability . We here evaluate LPJmL4 by comparing the simulated and observed yield variability of the 10 top-producing countries of the respective crop . We refrain from comparing to individual sentinel sites, but refer to the evaluation of LPJmL crop simulations at global, national, and grid cell scale in the global gridded crop model evaluation framework . As in , we aggregate simulated grid-cell-level yield time series to average national yield time series using the MIRCA2000 dataset for spatial aggregation and removing trends in observations and simulations with a moving-window average (see , for details). The productivity of biomass plantations is evaluated with data from experimental sites for Miscanthus, switchgrass, poplar, willow, and eucalyptus production using the data collection of . Data on biomass productivity typically report a data range. These are site-specific management differences and reflect the diverse drivers of reported productivity, such as variation in plant species, fertilizer use and irrigation management, crop spacing, or sapling size. We average the minimum and maximum values to derive the mean productivity per site.

Sowing dates

To evaluate the accuracy of the simulated rain-fed sowing dates, we use the global dataset of growing areas and growing periods, MIRCA2000 , at a spatial resolution of 0.5${}^{\circ }$ and a temporal resolution of 1 month, as proposed by . Monthly data in MIRCA2000 were converted to daily data by assuming that the growing period starts on the first day of the month following . MIRCA2000 reports several growing periods in a year for some administrative units for the crops wheat, rapeseed, rice, cassava, and maize. For comparison we select the best corresponding growing period so that a close agreement indicates that simulated sowing dates are reasonable, but not necessarily the most frequently chosen by farmers. We do not compare simulated sowing dates for sugar cane (see Fig. S94 in the Supplement) to observed sowing dates, as MIRCA2000 assumes it is grown all year-round as a perennial crop.

Evaluation metrics

We employ Taylor diagrams to compare the correlation, differences in standard deviation, and the centred root mean squared error (CRMS) between simulated and observed carbon and water fluxes at FLUXNET sites and at gauge stations from ArcticNET and UNH/GRDC. The standard deviations of the reference datasets have been normalized to 1.0 so that multiple sites can be displayed in one figure.

For global gridded reference datasets, such as for carbon stocks, we show spatial patterns in maps and aggregations as latitudinal means and quantify overall differences as a spatial correlation analysis over all grid cells (see Table ). As suggested by we use the normalized mean square error (NMSE) to describe differences between model simulation and reference datasets. The NMSE is zero for perfect agreement, 1.0 if the model is as good as using the data mean as a predictor, and larger 1.0 if the model performs less well than that. The squared error term puts stronger emphasis on large deviations between simulations and observations and is thus stricter than the normalized mean error (see Table  for equations). also suggest using the normalized mean error (NME) as a more robust metric than NMSE. NME is based on absolute residuals (NMSE on squared residuals) and thus is especially better suited for variables that can have very large values and residuals. Additionally, we use the Manhattan metric (MM) proposed by for evaluation of vegetation cover. Values for MM less than 1 reflect the fact that the model performs better than the mean value. Additionally, we show the random model, which was generated by bootstrap resampling of the observations as proposed by Table 4. The random model was used for the evaluation of vegetation distribution. Table  gives an overview of variables evaluated at the local scale and the measures that were used for the evaluation of time series for crop yields. We employ a simple time series correlation analysis after removing trends with a moving-window detrending method. For comparison with point measurements, we extract the time series from corresponding 0.5${}^{\circ }$ grid cells. These simulated time series may differ in terms of weather and soil conditions from the actual site as the simulations are based on gridded global dataset inputs. The time period is given by the respective measurements, which differ for each observation point.

Evaluation metrics used in this study.

Metric	Equation	Reference
NMSE	$\mathrm{NMSE}=\frac{{\sum }_{i=\mathrm{1}}^N(y_i-x_i{)}^{\mathrm{2}}}{{\sum }_{i=\mathrm{1}}^N(x_i-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}}$
NME	$\mathrm{NME}=\frac{{\sum }_{i=\mathrm{1}}^N\left|y_i,-,x_i\right|}{{\sum }_{i=\mathrm{1}}^N\left|x_i,-,\stackrel{\mathrm{‾}}{x}\right|}$
ME	$\mathrm{ME}=\frac{{\sum }_{i=\mathrm{1}}^N\left|y_i,-,x_i\right|\cdot A_i}{{\sum }_{i=\mathrm{1}}^N{A}_i}$
W	$\mathrm{W}=\mathrm{1}-\frac{{\sum }_{i=\mathrm{1}}^N{\left(y_i,-,x_i\right)}^{\mathrm{2}}\cdot A_i}{{\sum }_{i=\mathrm{1}}^N{\left(\left|y_i,-,\stackrel{\mathrm{‾}}{x}\right|,+,\left|x_i,-,\stackrel{\mathrm{‾}}{x}\right|\right)}^{\mathrm{2}}\cdot A_i}$
MM	$\mathrm{MM}=\frac{{\sum }_{i=\mathrm{1}}^N\left|q_{i,j},-,p_{i,j}\right|}{N}$

Note: $y_i$ is the simulated and $x_i$ the observed value in grid cell $i$, $\stackrel{\mathrm{‾}}{x}$ the mean observed value, $A_i$ the area weight in grid cell $i$, and $N$ the number of grid cells or sites; $q_{i,j}$ is the simulated and $p_{i,j}$ is the observed fraction of item $j$ in grid cell $i$. Normalized mean square error – NMSE, normalized mean error – NME, ME – mean absolute error, W – Willmott coefficient of agreement, MM – Manhattan metric.

Overview of variables and measures used for the evaluation of LPJmL4 local scale.

	Measure				Reference
		Standard		Reference
Variable	CRMSE	deviation	Correlation	to figures	Data	Citation
CO${}_{\mathrm{2}}$			x	Figs. ,	Atmospheric transport${}^{\mathrm{1}}$
NEE	x	x	x	Fig.	FLUXNET${}^{\mathrm{2}}$
ET	x	x	x	Fig.	FLUXNET${}^{\mathrm{2}}$
NPP			x	Fig. d
GPP			x	Fig. c
BIOMASS			x	Fig. a, b
DISCHARGE	x	x	x	Figs. ,	ArcticNET${}^{\mathrm{3}}$,
				S19–S66	UNH/GRDC${}^{\mathrm{4}}$

Centred root mean square error (CRMSE). ${}^{\mathrm{1}}$ http://pubman.mpdl.mpg.de/pubman/faces/viewItemOverviewPage.jsp?itemId=escidoc:1691952; ${}^{\mathrm{2}}$ http://fluxnet.fluxdata.org/data/la-thuile-dataset/; ${}^{\mathrm{3}}$ http://www.r-arcticnet.sr.unh.edu/v4.0/index.html; ${}^{\mathrm{4}}$ http://www.grdc.sr.unh.edu/index.html.

To envisage the degree of agreement between simulated (LPJmL4) and observed (MIRCA2000) sowing dates, we follow and compute two different metrics: the Willmott coefficient of agreement (W) and the mean absolute error (ME), both weighted by the crop-specific cultivated area according to . For an overview of all metrics used, see Table .

Results and discussion

In the following we compare the standard version LPJmL4, which refers to the experiment LPJmL4-GSI-GlobFIRM. In the case of the other experiments we refer to the names defined in Sect. .

Vegetation cover

LPJmL4 reproduces the observed vegetation distribution better than the random model (Table ). LPJmL4 can best reproduce the distinction between bare soil and vegetated areas (MM $=$ 0.22) and between tree-covered areas and areas without trees (MM $=$ 0.31), but with considerably better scores than the random model (MM $=$ 0.56 and 0.54, respectively). Moreover, LPJmL4 simulation results reach the lowest MM scores for the distinction of evergreen vs. deciduous trees (MM $=$ 0.52) and for the distribution and composition of life forms (trees vs. herbaceous vs. bare soil; MM $=$ 0.45); these are substantially better than the random model (MM $=$ 0.87 and 0.88, respectively). The largest improvement in LPJmL4 simulations over the random model is found for the patterns of broadleaved vs. needle-leaved trees (MM $=$ 0.37 for LPJmL4 vs. 0.94 for the random model; see Table ).

Comparison metric scores for LPJmL4 simulations against observations of fractional vegetation cover data from International Satellite Land-Surface Climatology Project (ISLSCP) II vegetation continuous field (VCF) .

Vegetation cover	Manhattan metric (MM)
	LPJmL4	Random model${}^{*}$
Life forms	0.45	0.88
Tree vs. non-tree	0.31	0.54
Herb vs. non-herb	0.42	0.66
Bare vs. covered ground	0.22	0.56
Evergreen vs. deciduous	0.52	0.87
Broadleaf vs. needle-leaf	0.37	0.94

MM suggested by , ${}^{*}$ values taken from Table 4.

Comparison of the atmospheric CO${}_{\mathrm{2}}$ concentrations at Point Barrow (BRW; panel a) and Mauna Loa (MLO; panel b) for the different LPJmL4 experiments.

[Figure omitted. See PDF]

Atmospheric CO${}_{\mathrm{2}}$ concentration and NEE

Comparison of the atmospheric CO${}_{\mathrm{2}}$ concentration at Mauna Loa (MLO) and Point Barrow (BRW) simulated in the different LPJmL4 experiments. (a, b) Seasonal cycle, (c) trend of the seasonal amplitude, and slopes are given for the different LPJmL4 experiments.

[Figure omitted. See PDF]

Comparison of simulated NBP to atmospheric CO${}_{\mathrm{2}}$ concentration at MLO and BRW

LPJmL4 reproduces the observed long-term and seasonal dynamics of atmospheric CO${}_{\mathrm{2}}$ well (Figs.  and ). The long-term trend of atmospheric CO${}_{\mathrm{2}}$ is reproduced well in all the different model set-ups (Fig. ), except for the set-up with natural vegetation only (LPJmL4-GSI-GlobFIRM-PNV). The experiment with all processes included (LPJmL4-GSI-GlobFIRM) gives the best correlation and trend reproduction, which suggests that an integral representation of the LPJmL4 features is required to match observations best. Next to land use dynamics, the inclusion of permafrost dynamics has the strongest effects on the simulated trend (LPJmL4-NOGSI-NOPERRM-GlobFIRM vs. LPJmL4-NOGSI-GlobFIRM). The use of the process-based fire model SPITFIRE leads to a small overestimation of the trend in atmospheric CO${}_{\mathrm{2}}$ concentrations compared to the other model set-ups, especially at MLO. Seasonal variations in atmospheric CO${}_{\mathrm{2}}$ can be reproduced well by LPJmL4, especially by the standard set-up (LPJmL4-GSI-GlobFIRM; Fig. ). The simulation of seasonal variations in atmospheric CO${}_{\mathrm{2}}$ content are especially improved by the GSI phenology scheme (LPJmL4-NOGSI-GlobFIRM vs. LPJmL4-GSI-GlobFIRM; Fig. a, b). All model set-ups (except LPJmL4-GSI-SPITFIRE) can reproduce the observed strong significant increase in the seasonal CO${}_{\mathrm{2}}$ amplitude at BRW and the weak (and insignificant) increase at MLO (Fig. c). These results are in agreement with a previous evaluation of simulated seasonal CO${}_{\mathrm{2}}$ changes in LPJmL .

Further analysis shows that the standard set-up (LPJmL4-GSI-GlobFIRM) can best produce the mean seasonal cycle in MLO, whereas the version that omits land use (LPJmL4-GSI-GlobFIRM-PNV) performs slightly better than this in BRW (Fig. ). The standard set-up (LPJmL4-GSI-GlobFIRM) can also best reproduce the increase in the seasonal amplitude at BRW, whereas it is the only set-up that produces a statistically significant but still very small increase in the seasonal amplitude at MLO where observations also do not show a statistically significant increase.

Comparison of simulated NEE to eddy flux measurements

Net ecosystem exchange rate measured at eddy flux towers: . Available online at FLUXNET (http://fluxnet.fluxdata.org/data/la-thuile-dataset/). Sites (colours) are ordered from north to south.

[Figure omitted. See PDF]

Overview of variables evaluating LPJmL4 showing measures and references at the global scale.

	Measure					Reference
			Spatial	Temporal	Visual
Variable	NME	NMSE	correlation	correlation	comparison	Data	Citation
GPP – Av	0.20	0.13	0.87		Figs. ,	GPP${}^{\mathrm{1}}$
					S68
$R_{\text{e}}$ – Av	0.67	0.55	0.67		Figs. ,
					S70
SoilC – Av	0.48	0.75	0.29		Fig. S67	Soil carbon stocks${}^{\mathrm{1}}$
VegC – Av	0.33	0.36	0.84		Fig. S69a	Total biomass${}^{\mathrm{1}}$
					Fig. S69b	AGB
FAPAR – I-aMv	0.17	0.13	0.63	Fig. a		MODIS FAPAR${}^{\mathrm{2}}$
FAPAR – I-aMv	0.18	0.15	0.59	Fig. b		GIMMS3g FAPAR${}^{\mathrm{3}}$
FAPAR – I-aMv	0.21	0.20	0.69	Fig. c		VGT2 FAPAR${}^{\mathrm{4}}$
ET	1E-6	0.07	0.84		Fig. S71	Latent heat flux${}^{\mathrm{1}}$
fBA					Fig. S72		GFED4 & CCI Fire (4.1)
Albedo					Fig. S72	MODIS C5
Discharge						ArcticNET${}^{\mathrm{5}}$,
Ov	0.42	0.24		$R^{\mathrm{2}}=\mathrm{0.90}$		UNH/GRDC${}^{\mathrm{6}}$
Mav	0.36	0.19		$R^{\mathrm{2}}=\mathrm{0.92}$
I-av	0.24	0.06		$R^{\mathrm{2}}=\mathrm{0.97}$

Normalized mean error (NME) and normalized mean square error (NMSE) as suggested by ; Av – annual average; I-aMv – inter-annual monthly variability; overall variability – Ov; monthly average variability – Mav; inter-annual variability – I-av; vegetation carbon – VegC; aboveground biomass – AGB; soil carbon – SoilC; fBA – fractional burnt area.${}^{\mathrm{1}}$ https://www.bgc-jena.mpg.de/geodb/BGI/Home; ${}^{\mathrm{2}}$ https://lpdaac.usgs.gov/dataset_discovery/modis/modis_products_table/mod15a2; ${}^{\mathrm{3}}$ http://cliveg.bu.edu/modismisr/lai3g-fpar3g.html; ${}^{\mathrm{4}}$ http://cordis.europa.eu/result/rcn/140496_en.html; ${}^{\mathrm{5}}$ http://www.r-arcticnet.sr.unh.edu/v4.0/index.html; ${}^{\mathrm{6}}$ http://www.grdc.sr.unh.edu/index.html.

We evaluate the model performance of simulated NEE from LPJmL4 for temporal and spatial variation in NEE data from eddy flux measurements using Taylor diagrams . Stations are sorted from north to south (see Fig. ) for all NEE measurements available for $>$ 3 years. The model is able to reproduce the mid-latitudes best (represented by yellow over green to light blue colours), with correlation coefficients mostly between 0.4 and 0.9 and standard deviations often within $\pm$30 $\mathrm{\%}$ of the reference data. The northernmost regions are reproduced well at some flux towers, but often with higher standard deviation than in the flux tower data, which means that the simulated time series are largely in phase but more variable than the observations. In contrast, the evaluation is comparatively poor for tropical regions, especially the station at Santarém with strong negative correlations ($r<-\mathrm{0.6}$) but realistic standard deviations. For this site, however, have already pointed out that the eddy flux measurements show the opposite sign compared to tree growth observations and model predictions, which is also the case for LPJmL4. We stress that this evaluation is done for a standard LPJmL4 run and standard input (the LPJmL4-GSI-GlobFIRM as described in ); i.e. we did not calibrate the model to site-specific conditions and also drive the model with gridded input data rather than the observed soil and weather data at individual stations. More detail for comparisons with eddy flux tower measurements for individual locations is supplied in the Supplement (see Figs. S1–S7). Additionally, we have simulated NEE by conducting simulations with station-specific meteorological observations (see Fig. S17). The results are similar to simulations driven by global climate data.

Vegetation and soil carbon stocks and vegetation productivity

Soil carbon and vegetation carbon stocks

The spatial correlation between simulated and observation-based estimates of SOC by is weak ($r=\mathrm{0.29}$; Table ) with disagreements in the subtropics where LPJmL4 simulations substantially underestimate soil carbon stocks, whereas LPJmL4 reports much higher soil carbon in the high northern latitudes ($>$ 50${}^{\circ }$ N) and lower values for the tropical and temperate zone compared to (see Fig. S67). Other estimates by show much higher carbon content for the permafrost-affected areas than the dataset of . We thus assume that the disagreement between simulations and the data may also result from an underestimation of carbon stocks in the data. However, the estimation of global soil carbon is less in LPJmL4 (1869 $\mathrm{Pg}\mathrm{C}$) than estimated by (2352 $\pm$ 400 $\mathrm{Pg}\mathrm{C}$).

Evaluation of vegetation carbon (a), aboveground biomass (b), GPP (c), and NPP (d). Observed data are provided by . Bars give the minimum and maximum of the estimation within one 0.5${}^{\circ }$ cell simulated by LPJmL4.

[Figure omitted. See PDF]

The comparison of simulated and observation-based assessments of vegetation carbon show a good spatial correlation ($r=\mathrm{0.84}$; Table ). Globally, estimate slightly lower biomass (445 $\pm$ 8 $\mathrm{Pg}\mathrm{C}$) as simulated by LPJmL4 (507 $\mathrm{Pg}\mathrm{C}$). The spatial patterns of vegetation carbon stocks are shown in Fig. S69a for simulations and the data product of . While the broad geographical patterns are in overall agreement with the evaluation data, the absolute values differ in some regions. Specifically, LPJmL4 simulates much higher biomass (see the latitudinal pattern in Fig. S69) for the tropics and lower biomass between 20 and 40${}^{\circ }$ in the Northern and Southern Hemisphere where show higher values compared to LPJmL4. This is probably due to an overestimation of vegetation carbon in agricultural regions by , as show similar aboveground biomass estimates there (see Fig. S69b). The subtropical region where biomass carbon is underestimated also corresponds to the region where LPJmL4 simulations underestimate soil carbon stocks compared to . Also, the comparison of aboveground biomass estimates with the dataset of shows a similar spatial pattern of overestimation of vegetation biomass with too-high values in boreal and tropical areas. The comparison is complicated by uncertainties in the estimation of belowground biomass and the assumed distribution between aboveground and belowground biomass in LPJmL4 simulations, in which LPJmL4 assumes that belowground biomass consists of all fine root biomass and one-third of all sapwood biomass. The simulation experiments without permafrost dynamics (LPJmL4-NOGSI-NOPERM-GlobFIRM) show a high overestimation of biomass in the high latitudes. Similarly, the inclusion of the GSI phenology substantially reduces the biomass overestimation in comparison to and , which is consistent with the finding of . The consideration of human land use in the simulations improves carbon stock simulations in the temperate zones (Fig. S69). This clearly demonstrates the importance of permafrost, human land use, and the GSI phenology for the simulation of the terrestrial carbon cycle, even though the remaining discrepancies warrant further model improvement.

The maps (a) show the spatial pattern of gross primary production (GPP; $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{a}}^{-\mathrm{1}}$) distribution from the standard LPJmL4 simulation against the MTE data . The graph in (b) shows the latitudinal pattern of GPP distribution simulated by the different versions of LPJmL4 against data from .

[Figure omitted. See PDF]

Figure a and b compare site data estimation with the representative LPJmL4 grid cell estimation with an uncertainty range which comes from the different measurements within one 0.5${}^{\circ }$ grid cell. Both vegetation and aboveground carbon are slightly overestimated in some cases but also strongly underestimated in others. As LPJmL4 calculates a representative mean value of a 0.5${}^{\circ }$ grid cell for all benchmarks, the simulated values should match the mean values. However, it can be assumed that measurements are not evenly distributed through the age classes within one grid cell or forest, and it remains unclear how representative the measurements are for a 0.5${}^{\circ }$ grid cell area.

Gross and net primary production (GPP and NPP)

The global estimation of 123.7 $\mathrm{Pg}\mathrm{C}{\mathrm{a}}^{-\mathrm{1}}$ GPP from LPJmL4 (see Fig. ) matches the estimates from and of 123 $\pm$ 8 and 119 $\pm$ 6 $\mathrm{Pg}\mathrm{C}{\mathrm{a}}^{-\mathrm{1}}$, respectively, for the years 1982–2005, whereas the highest divergence can be observed in the tropics where LPJmL4 estimates much lower values despite the higher biomass estimations (see Sect. ). LPJmL4 simulated higher GPP for the temperate and boreal zones than reported by . The different model experiments show similar patterns except for LPJmL4-GSI-GlobFIRM-PNV, which shows lower GPP in the Mediterranean (see Fig. ). estimate global NPP at 54 $\pm$ 10 $\mathrm{Pg}\mathrm{C}{\mathrm{a}}^{-\mathrm{1}}$ and LPJmL4 at 57 $\mathrm{Pg}\mathrm{C}{\mathrm{a}}^{-\mathrm{1}}$ for the mean of the years 1982–2011.

The site data comparison to shows a good agreement between site measurements and simulated GPP (see Fig. c) and NPP (see Fig. d). The overestimation of simulated biomass and the good agreement of NPP and GPP leads to the conclusion that LPJmL4 underestimates mortality. This warrants further investigation of why LPJmL4 seems to overestimate global GPP but shows good agreement with site data. The comparison of LPJmL4 against MTE data on the local scale for the same points as given by shows a good agreement, especially if outliers are excluded (Fig. S68b, c). Figure S68a compares plot data against the global data.

Ecosystem respiration ($R_{\text{e}}$)

Comparison of satellite-derived ecosystem respiration with that simulated by LPJmL4 reveals similar spatial patterns (Figs.  and S70). However, LPJmL4 shows higher temperature sensitivities (Fig. a) and consistently simulates higher $R_{\text{e}}$ values in high-latitude and subtropical regions (Fig. S70). Since satellite-derived ecosystem respiration is calibrated for FLUXNET data and hence exhibits marginal cross-latitude bias, the discrepancies to LPJmL4 are likely associated either with LPJmL4 parameterization or with systematic errors in the FLUXNET processing technique. Additional details and figures are presented in .

Ecosystem respiration ($R_{\text{e}}$) evaluation of standard LPJmL4 simulations with satellite-derived estimations from . Annual $R_{\text{e}}$ sums for all pixels from the displayed extent in Fig. S70 are compared and separated by climate type (a)–(c). Dashed lines indicate a polynomial bias curve. Chart symbols are separated for forest (FOR) and grassland–cropland (GRA–CRO) land cover classes.

[Figure omitted. See PDF]

Water fluxes

Evapotranspiration

The spatial distribution of evapotranspiration in LPJmL4 shows a very similar pattern to that estimated by (Table , Fig. S71). It indicates a general underestimation of ET, especially in the tropics and subtropics, but in most cases within the uncertainty range. This is consistent with the underestimation of GPP in the tropics (Fig. ), but not with the general overestimation of vegetation biomass (Fig. S69). The different experiments show nearly no effects on the simulated evapotranspiration. At site level, the evapotranspiration fluxes show a good agreement with eddy flux tower measurements (Fig. ). LPJmL4 shows good performance in most regions, with correlation coefficients often larger than 0.6. The northern and temperate stations (red to light blue symbols) show especially high correlation with low CRMS. Simulations of tropical and subtropical ET (dark blue to purple symbols) show weak or even negative correlations coupled with a high CRMS for some stations. We also provide more detailed time series analyses for the evapotranspiration fluxes of individual sites in the Supplement (Figs. S8–S16).

Evaporation rate measured at eddy flux towers: . Available online at FLUXNET (http://fluxnet.fluxdata.org/data/la-thuile-dataset/). Site locations are ordered from north to south.

[Figure omitted. See PDF]

Comparison of simulated discharge with 287 gauges provided by ArcticNET (http://www.r-arcticnet.sr.unh.edu/v4.0/index.html) and UNH/GRDC (http://www.grdc.sr.unh.edu/index.html). Stations with basin area $\ge$ 10 000 ${\mathrm{km}}^{\mathrm{2}}$ are taken into account. Gauges are ordered from north to south (see legend colour).

[Figure omitted. See PDF]

River discharge stations evaluation

Discharge simulated by earlier LPJmL versions was previously evaluated in several studies, also in comparison with other global hydrological and land surface models . River discharge was evaluated for major catchments globally, also accounting for the effects of different precipitation datasets and regionally for the Amazon basin and the Ganges . Figure  shows the comparison of simulated LPJmL4 and observed river discharge values for all gauges with a basin area $\ge$ 10 000 ${\mathrm{km}}^{\mathrm{2}}$. Here, the most northern (blue) and also the most southern (purple) gauges show good agreement, but overall the picture is mixed with respect to correlation coefficients and standard deviation. For further insights, we provide comparisons for all considered gauges in the Supplement (Figs. S19–S66). For many gauges, the simulated seasonal timing of river discharge (peaks) has improved (see Figs. S19–S22) compared to the previous model evaluation of river discharge , which is mainly a result of the newly implemented GSI phenology scheme . The discharge spring peaks in permafrost areas are especially affected by this improvement. At many gauges, LPJmL4 can reproduce the variability for the whole time series and specially the seasonality, with a high $R^{\mathrm{2}}$ and NME and NMSE, which implies a better performance than the mean model. The dynamics at gauges in the temperate zone (Figs. S49–S50, S61) are not well reproduced in the simulations, and the NME and NMSE also show high values in contrast to gauges in the subtropics and tropics (Figs. S64–S66), which typically show high $R^{\mathrm{2}}$ and low NME and NMSE.

The evaluation at the global aggregation (computed for all stations and then averaged) shows very high agreement between observed and modelled discharge (see Table ). Both the explained variance ($R^{\mathrm{2}}$) and the NME–NMSE contribute to the good performance of the simulated discharge. The constant flow velocity in all rivers, as assumed in LPJmL4 simulations, could be varied by river for further model improvement, especially for the timing in flat areas where wetland dynamics may play an important role.

Irrigation withdrawal and consumption

Global estimates of irrigation water withdrawal ($W_{\text{d}}$: 2545 ${\mathrm{km}}^{\mathrm{3}}{\mathrm{a}}^{-\mathrm{1}}$) and consumption ($W_{\text{c}}$: 1292 ${\mathrm{km}}^{\mathrm{3}}{\mathrm{a}}^{-\mathrm{1}}$) agree well with previous studies. Reported $W_{\text{d}}$ values for the period 1998–2012 are 2722 ${\mathrm{km}}^{\mathrm{3}}{\mathrm{a}}^{-\mathrm{1}}$ , and modelling results range from 2217 to 3185 ${\mathrm{km}}^{\mathrm{3}}{\mathrm{a}}^{-\mathrm{1}}$ . $W_{\text{c}}$ estimations range between 927 and 1530 ${\mathrm{km}}^{\mathrm{3}}{\mathrm{a}}^{-\mathrm{1}}$ . find that 1179 ${\mathrm{km}}^{\mathrm{3}}{\mathrm{a}}^{-\mathrm{1}}$ (1098 ${\mathrm{km}}^{\mathrm{3}}{\mathrm{a}}^{-\mathrm{1}}$ in ) relates to surface water with an additional 257 ${\mathrm{km}}^{\mathrm{3}}{\mathrm{a}}^{-\mathrm{1}}$ from groundwater resources. LPJmL4 does not account for fossil groundwater extraction nor desalination. However, previous studies show that 80 % of groundwater withdrawals are recharged by return flows . It is thus plausible that studies accounting for (fossil) groundwater reach $W_{\text{d}}$ estimates somewhat higher than in LPJmL4. Naturally, irrigation water estimates are associated with uncertainties in the precipitation input employed . A representation of multiple cropping systems in LPJmL4 and the corresponding growing seasons could also help to improve water withdrawal and consumption estimates and eventually river discharge, especially in tropical areas.

Simulated irrigation efficiencies are difficult to compare with observations due to inhomogeneous definitions and field measurement problems. Yet, in Table S1 in the Supplement we relate our results to comparable literature. Our simulations meet the indicative estimates of at the global level. provide another independent estimate of field efficiency with global average values of 42, 78, and 89 $\mathrm{\%}$ for the three irrigation types, respectively. Our estimates agree well with these numbers globally and regionally, even though there are some regional patterns that are not represented in our results. , for instance, find lower surface irrigation efficiencies in the Middle East, North Africa (MENA), and sub-Saharan Africa (SSA). We simulate above-average efficiencies in MENA and particularly low ones in South Asia, which are both supported by and . Overall, the evaluation of the irrigation model in LPJmL4 demonstrates that it is well in line with reported patterns, and yet it comes with much more detailed depths with respect to process representation and spatio-temporal resolution than these.

Permafrost distribution and active layer thickness

The current permafrost distribution and the active layer thickness (Fig. ) is well represented by the LPJmL4 model compared to independent studies . LPJmL4 is able to reproduce the distribution of permafrost and the measured active layer thickness in most grid cells. The continuous permafrost zone is characterized by a thawing depth equal to or less than 1 m in LPJmL4, while the model simulates for sporadic permafrost and isolated patches a thawing depth of more than 3 m. The spatial distribution of greater thaw depth from north to south is simulated well by the model. CALM station data show a similar thawing depth as simulated by LPJmL4 (Fig. b), but CALM station data also indicate that thawing depth can be different for the same grid cell, as other processes (e.g. exposition) not represented by LPJmL4 can play an important role.

Observed and simulated permafrost distribution and active layer thickness. (a) Contemporary permafrost extent according to the IPA Circum-Arctic Map of Permafrost (${}^{*\mathrm{1}}$). (b) LPJmL4-simulated active layer thickness compared to the ${}^{*\mathrm{2}}$CALM station data means both for the observation time 1991–2009 (https://www2.gwu.edu/~calm/; ). The colour scheme used at the bottom is the same for simulated thaw depth and Circumpolar Active Layer Monitoring (CALM) data.

[Figure omitted. See PDF]

Fire

Burnt area

Simulated fractional area burnt is largest in the seasonal dry tropics and temperate regions in all model versions and smallest in cold or wet environments (Fig. S72). However, maximum fractional burnt area does not exceed 0.0625 in tropical and subtropical savanna and shrubland areas when the GlobFIRM model is applied. It is comparable to GFED4 and CCI estimates only in South America, while in other tropical regions GFED4 and CCI report fractional burnt area between 0.125 and 0.75 (Fig. S72). In these regions, the fractional burnt area simulated by the SPITFIRE model is overestimated with values between 0.25 and 1, specifically in Southern Hemispheric Africa and northern Australia. SPITFIRE is very sensitive to vegetation, and thus fuel composition in areas with homogeneous C${}_{\mathrm{4}}$ grasslands can lead to an overestimation of simulated area burnt, which is specifically the case for seasonally dry South America and the Indian subcontinent. LPJmL4-GSI-SPITFIRE captures the distribution of fractional burnt area much better than LPJmL4-GSI-GlobFIRM, which is too homogeneous in its response. In contrast, LPJmL4-GSI-SPITFIRE better captures the very small fractions reported for the wet tropical forests, which is better comparable to GFED4. Here, the approach of simulating fire risk based on the climatic fire danger index instead of deriving a fire probability from the topsoil soil moisture is of great advantage in these regions. While LPJmL4-GSI-GlobFIRM simulates a relatively homogeneous spatial distribution of fractional burnt area in temperate and boreal forest regions, LPJmL4-GSI-SPITFIRE underestimates fractional burnt area in these biomes. LPJmL4-GSI-GlobFIRM underestimates fractional burnt area in the temperate steppe regions, whereas LPJmL4-GSI-SPITFIRE manages to spatially capture the burning conditions in these biomes, even though the total amount is overestimated. The phenology module in LPJmL4 has no effect on the fractional burnt area simulated by LPJmL4-GSI-GlobFIRM, whereas including permafrost increases burnt area in the circum-boreal region, specifically in Siberia, even though the spatial effect is too homogeneous.

Fire effects on biomass and vegetation distribution

Both fire model approaches simulate a comparable latitudinal distribution of biomass starting from the wet tropics towards dry and colder areas in the north and south. Both model versions simulate comparable values in the wet tropics around the Equator and capture the gradient to seasonal dry tropics in the north (until 10${}^{\circ }$ N) and south (until 20${}^{\circ }$ S). The overestimation of burnt area in tropical savannas around 20${}^{\circ }$ N in LPJmL4-GSI-SPITFIRE leads to an underestimation in simulated biomass compared to the other LPJmL4 experiments. The consideration of permafrost and fire dynamics is required to reproduce observed vegetation biomass values in boreal regions.

Global biomass burning

The modelling errors in fractional area burnt compensate in different ways in each fire model. SPITFIRE simulates global biomass burning values of 2.7 $\mathrm{Pg}\mathrm{C}{\mathrm{a}}^{-\mathrm{1}}$ on average between 1996 and 2005, which is comparable to the 2.33 $\mathrm{Pg}\mathrm{C}{\mathrm{a}}^{-\mathrm{1}}$ suggested by . Here, overestimations of burnt area in tropical savannas and underestimations in boreal forests compensate for each other. GlobFIRM simulates more fires in boreal regions that are less spatially pronounced than in GFED4, but underestimates fractional burnt area in the subtropics and tropics. GlobFIRM therefore estimates global biomass burning by 2.8 $\mathrm{Pg}\mathrm{C}{\mathrm{a}}^{-\mathrm{1}}$, which is similar to SPITFIRE.

Fraction of absorbed photosynthetically active radiation (FAPAR) and albedo

Evaluations against multiple satellite datasets of FAPAR have already shown that LPJmL-GSI can reproduce the seasonality of FAPAR and the inter-annual variability and trends well at the start and end of the growing season within observational uncertainties . LPJmL4 shows a high spatial correlation with correlation coefficients between 0.6 and 0.71 for PEAK-FAPAR. It shows also a good agreement with the temporal variations (Fig. a–c). Large parts of the wet tropics display a negative correlation between simulated and observed FAPAR, which may explain the phase offset in the dynamics of NEE at the station Santarém. However, in these regions the differences between datasets are also large, which is caused by the limitations of optical satellite observations in regions with permanent cloud cover .

Evaluation of FAPAR for different data sources: MODIS (a), GIMMS (b), and VGT2 (c).

[Figure omitted. See PDF]

LPJmL4 reproduces the global patterns of annual peak FAPAR (Fig. ) well. In northern latitudes and in the tropics, LPJmL4 is within the range of the FAPAR datasets. However, LPJmL4 overestimates peak FAPAR, especially in middle and low latitudes, which originates from an overestimation of FAPAR in semi-arid regions. LPJmL4 reproduces the temporal dynamic of FAPAR well in most climate regions with very high correlations between simulated and observed FAPAR in temperate and boreal climates (climate regions Cf and D${}^{*}$) and with medium to high correlations in semi-arid climate regions (e.g. Am, As, Aw, Bsh, Bsk, Cs in Fig. S73). LPJmL4 and the observational datasets show low correlations in wet tropics (Af) and in winter-dry temperate climates (Cw).

FAPAR mean annual peak comparison with three different remote sensing products.

[Figure omitted. See PDF]

LPJmL4 overestimates albedo in all regions (Fig. S74). The temporal dynamic of snow-free albedo was reproduced well in cold steppes (climate region BSk) and in boreal regions (climate regions D${}^{*}$). The correlation between simulated and observed albedo is poor in tropical semi-arid and temperate climates (e.g. As, Aw, Cs, Cf). This is likely caused by soil-moisture-induced changes in soil and background albedo, which has a great effect on soil reflectance outside the vegetation season. Such changes are not considered in LPJmL4.

Agriculture

Crop yields variability

The evaluation of simulated crop growth and yield can be assessed at individual sites if the model is used as a point model as in different model intercomparison simulations in which reference data are available for end-of-season properties (most importantly, crop yield) and within-season dynamics (e.g. development of leaf area index, LAI). The crop yield simulations of LPJmL were evaluated in the framework of the Agricultural Model Intercomparison and Improvement Project (AgMIP) for wheat, maize, rice, and soybean by . They find that the performance of LPJmL is similar to that of the other gridded crop models in that model ensemble ($n=\mathrm{14}$). We supplement the model evaluation with time series correlation analyses for the 10 top-producing countries for all crops implemented in LPJmL4 . Results are portrayed in Fig. , except for field peas for which no spatial data on crop-specific harvested areas exist for aggregation to national yield time series . As national yield levels are roughly calibrated in standard LPJmL simulations , a comparison of the mean bias does not provide insights on model performance. As management intensity is assumed to be static in the simulations (Sect. ), yield trends cannot be reproduced so that simulated and reported national yield time series have been detrended with a running mean approach prior to comparison. For a more comprehensive evaluation of LPJmL's performance in yield simulations, see .

Evaluation of simulated yield variability for wheat (a) and maize (b) in comparison to FAO data (http://www.fao.org/faostat/en/#data/QC).

[Figure omitted. See PDF]

Map of simulated biomass yields by LPJmL4 from rain-fed herbaceous (a) and woody (b) BFTs (averages 1994–2009). Dots indicate the location of the experimental sites and measured yield, with colours scaled to map colours. Scatterplots compare observed and simulated yields in the respective grid cells. Model uncertainty is derived from simulations with and without irrigation. Observation uncertainty reflects dependencies on plantation management (adapted from ).

[Figure omitted. See PDF]

The agreement between simulated and observed yields is not only dependent on model performance, but also on the aggregation mask used , assumptions on management and model parameterization , soil parameters , and weather data inputs . LPJmL4 yield simulations are typically correlated with national yield statistics for some of the 10 top-producing countries for each crop, but only for one country in the case of cassava (Brazil) and sugar cane (China; Figs.  and S75–S83 for the other crops).

Biomass yield

For the purpose of this evaluation, irrigated and rain-fed biomass plants were simulated to grow globally wherever biophysical conditions allow for sustained growth. The averaged simulated yields for the 16-year period (1994–2009) were compared to reported biomass yields of switchgrass, Miscanthus, poplar, willow, and eucalyptus plantations on experimental test sites located in the respective grid cell (Fig. ). Simulated yields are mostly within the range of observations for Miscanthus, poplar, willow, and eucalyptus, but mostly overestimate switchgrass productivity. Management options for BFTs implemented in LPJmL4 are limited to irrigation management (rain-fed and fully irrigated) because plant species and plantation characteristics (e.g. sapling size and crop spacing) are parameterized as a constant scenario setting and were not varied here. The differences between rain-fed and irrigated biomass yield simulations are depicted as vertical error bars in Fig. . The range of rain-fed vs. fully irrigated biomass yields represents an approximation of management uncertainty because simulated yields depend strongly on water availability. Nevertheless, the simulated yield range is likely to represent optimal field management for rain-fed and irrigated plantations as nutrient limitations are not taken into account in these simulations.

Indices of agreement between simulated (LPJmL4) and observed (MIRCA2000) sowing dates.

	All cells			Precipitation seasonality			Temperature seasonality
Crop	W (–)	ME ($\mathrm{days}$)	$N$	W (–)	ME (days)	$N$ ($\mathrm{\%}$)	W (–)	ME ($\mathrm{days}$)	$N$ ($\mathrm{\%}$)
Wheat	0.87	44	13 962	0.86	40	15	0.87	44	85
Rice	0.90	25	4995	0.90	24	82	0.87	28	18
Maize	0.88	37	16 333	0.89	37	48	0.85	36	52
Millet	0.89	17	7851	0.92	16	63	0.89	31	37
Pulses	0.63	69	14 712	0.61	80	48	0.84	37	52
Sugar beet	0.37	19	2918	0.24			0.37	19	100
Cassava	0.93	51	6082	0.93	51	83	0.95	57	17
Sunflower	0.92	25	5876	0.87	45	22	0.93	22	78
Soybean	0.94	36	8259	0.94	35	31	0.92	36	69
Groundnut	0.77	34	5642	0.71	36	81	0.96	20	19
Rapeseed	0.86	49	5680	0.36	135	13	0.92	37	87
Wheat (excl. Russia)	0.94	30	11511	0.86	40	18	0.94	29	82

Mean absolute error (ME) and the Willmott coefficient of agreement (W).

Evaluation of sowing dates for wheat. From (a) to (c): simulated (LPJmL4) sowing date, observed (MIRCA2000) sowing date, and difference between simulated and observed sowing date. Green colours (red colours) in the difference map indicate that simulated sowing dates are too late (too early) compared to observations. White indicates crop area with less than 0.001 % of the grid cell area. Regions without seasonality are not shown.

[Figure omitted. See PDF]

Month of sowing

The average mean error (ME) for all crops globally is smaller than 2 months, with the exception of pulses (Table ). For wheat (excl. Russia), millet, rice, sunflower, and sugar beet, the agreement between the simulated and observed timing of sowing is higher, with a difference of about 1 month. The Willmott coefficients (W) are high, indicating good agreement between observations and simulations (W $>$ 0.85) for all crops except pulses, sugar beet, and groundnut. Both measures indicate closer agreement for pulses, groundnut, sunflower, and rapeseed in temperate regions . Poor agreement, with differences between simulated and observed sowing dates of more than 5 months, is found for maize and cassava in South-east Asia and China (for maize in East Africa), for wheat in Russia, for pulses in South-east Asia, India, West and East Africa, the south-east region of Brazil, and southern Australia, for groundnut in India and Indonesia, and for rapeseed in southern Australia and southern Europe (for wheat, Fig. ; for the other crops, Figs. S84–S93). Divergences are also substantial for crops growing in the southern part of the Democratic Republic of Congo, in South-east Asia, and in tropical climates.

There are several reasons for these disagreements between sowing dates simulated solely using climate data and the global crop calendar; please see for a more detailed discussion. Firstly the crop varieties in the crop calendar and those simulated here differ, i.e. spring and winter varieties of wheat and rapeseed in temperate regions (e.g. in Russia). Secondly, multiple cropping in tropical regions with high cropping intensity and complex cropping systems is not considered here. Thirdly, we use only one global temperature threshold for simulating sowing temperatures, which is known to vary between regions, and lastly, there are other uncertainties in our method of simulating sowing dates and in the global crop calendar we use for comparison. We also neglect important factors such as the availability of labour and machinery, social customs, markets and prizes, and the demand for certain agricultural products at certain times in the year.

The comparison to the global crop calendar, however, shows that close agreement between simulated and observed sowing dates can be achieved with purely climate-driven rules for large parts of the Earth for wheat, rice, maize, millet, soybean, and sunflower, as well as for pulses and groundnut in temperate regions. For about 75 % of the global cropping area the difference between simulated and observed sowing dates is 2 months; with the exception of cassava and rapeseed, 80 % of the crop area displays a difference of only 1 month, which is the minimum possible difference as the crop calendar reports monthly sowing dates.

Conclusions

This article provides a comprehensive evaluation of the now launched version 4.0 of the LPJmL DGVM that includes an operational representation of agriculture. Unique in its combination of features, the LPJmL4 model enables the simulation of carbon and water fluxes linked to the dynamics of both natural and agricultural vegetation in a single, internally consistent framework. We show that the model has great strength in reproducing carbon fluxes, especially for NBP on the global scale and NEE on the local scale. But we are also able to show that water fluxes match well with other estimates. Both carbon and water fluxes are the link to many ecosystem processes that the model represents and therefore are very important for the understanding of its interrelation. In the agriculture sector we conclude that in regions with a strong weather signal the model is able to match annual yield variability. Nevertheless, in highly managed countries yield variability is not well reproduced by the LPJmL4 model. This can be explained by the absence of a management module in the model. By following suggestions for objective intercomparative benchmarking systems of multiple models with dedicated software , the evaluation takes into account a number of performance metrics, diagnostic plots, and a broad range of fundamental model features. This work thus goes well beyond earlier evaluations of DGVMs (see ) and of model evaluations published for earlier versions LPJmL or its modules.

Pending major model improvements – anticipated as part of forthcoming LPJmL versions – are the incorporation of a scheme for calculating groundwater recharge and storage, the representation of nitrogen cycling for both natural and agricultural landscapes, consideration of ozone effects on plants and of soil degradation, representation of wetlands with associated methane emissions, the continuous refinement of crop parameterization including multi-cropping and other management forms, and possibly a revised implementation of soil moisture (following e.g. ) and stomatal conductance (following e.g. ). As such improvements are expected to have significant effects on plant production and carbon and water fluxes, thus influencing the overall model performance, any future LPJmL version will routinely be subjected to the evaluation protocol used here and, if applicable, tested against other standardized inter-model benchmarks (including participation in model intercomparisons with evaluation of single components such as in ). Such continued model maintenance and benchmarking shall also keep pace with recent developments in observational and experimental data, ideally supporting the identification of key uncertainties in model performance see. Besides identifying features for future model improvement, we demonstrate the adequate performance of the LPJmL4 DGVM in terms of the simulation of long-term averages and also the temporal dynamics across biogeochemical, hydrological, and agricultural processes. This unique capacity renders the LPJmL4 model suitable for process-based analyses of biosphere dynamics including assessments of multi-sectoral impacts of climate change or other anthropogenic Earth system interferences.

The model code of LPJmL4 is publicly available through PIK's gitlab server at https://gitlab.pik-potsdam.de/lpjml/LPJmL, and an exact version of the code described here is archived under 10.5880/pik.2018.002 and should be referenced as . The output data from the model simulations described here are available at the research data repository http://dataservices.gfz-potsdam.de/portal/ under 10.5880/pik.2017.009 and can be referenced as .

The Supplement related to this article is available online at https://doi.org/10.5194/gmd-11-1377-2018-supplement.

The authors declare that they have no conflict of interest.

Acknowledgements

This study was supported by the German Federal Ministry of Education and Research (BMBF) project “PalMod 2.3 Methankreislauf, Teilprojekt 2 Modellierung der Methanemissionen von Feucht- und Permafrostgebieten mit Hilfe von LPJmL” (FKZ 01LP1507C). Matthias Forkel was funded by the TU Wien Wissenschaftspreis 2015 awarded to Wouter Dorigo. This work used eddy covariance data acquired and shared by the FLUXNET community, including these networks: AmeriFlux, AfriFlux, AsiaFlux, CarboAfrica, CarboEuropeIP, CarboItaly, CarboMont, ChinaFlux, Fluxnet-Canada, GreenGrass, ICOS, KoFlux, LBA, NECC, OzFlux-TERN, TCOS-Siberia, and USCCC. The ERA-Interim reanalysis data are provided by ECMWF and processed by LSCE. The FLUXNET eddy covariance data processing and harmonization were carried out by the European Fluxes Database Cluster, AmeriFlux Management Project, and Fluxdata project of FLUXNET, with the support of CDIAC and ICOS Ecosystem Thematic Center, and the OzFlux, ChinaFlux, and AsiaFlux offices. We thank the coordinators of the Circumpolar Active Layer Monitoring (CALM) programme for providing thaw depth, the National Snow & Ice Data Center for providing the Circum-Arctic Map of Permafrost, and the R-ArcticNET for discharge data. Furthermore, we thank Jena-BGI data for providing GPP, latent heat flux, total biomass, and soil carbon data. We also thank the providers of the evaluation data of FAPAR (GIMMS3g FAPAR, VGT2 FAPAR, MODIS FAPAR) and the remote sensing data of GFED4 and CCI for evaluating fractional burnt area. MODIS C5 albedo time series data product was retrieved from the online data pool courtesy of the NASA Land Processes Distributed Active Archive Center (LP DAAC), USGS/Earth Resources Observation and Science (EROS) Center, Sioux Falls, South Dakota . We thank the Climatic Research Unit for providing global gridded temperature input, the Global Precipitation Climatology Centre for providing precipitation input, and the coordinators of ERA-Interim for providing shortwave downward radiation and net downward longwave radiation. Furthermore, we thank the coordinators of MIRCA2000 for providing land use input. Finally, we thank Kirsten Elger for her great support in archiving data and the LPJmL4 code and two anonymous reviewers for their helpful comments on earlier versions of the paper. Edited by: Julia Hargreaves Reviewed by: two anonymous referees