# Assessing the representation of the Australian carbon cycle in global vegetation models

Lina Teckentrup, Martin G de Kauwe, Andrew J Pitman, Daniel S Goll, Vanessa Haverd, Atul K Jain, Emilie Joetzjer, Etsushi Kato, Sebastian Lienert, Danica Lombardozzi, et al.

► To cite this version:

Lina Teckentrup, Martin G de Kauwe, Andrew J Pitman, Daniel S Goll, Vanessa Haverd, et al. Assessing the representation of the Australian carbon cycle in global vegetation models. Biogeosciences, 2021, 18 (20), pp.5639-5668. 10.5194/bg-18-5639-2021. hal-03415730

HAL Id: hal-03415730

https://hal.science/hal-03415730

Submitted on 5 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Assessing the representation of the Australian carbon cycle in global vegetation models

Lina Teckentrup<sup>1,2</sup>, Martin G. De Kauwe<sup>3,1,2</sup>, Andrew J. Pitman<sup>1,2</sup>, Daniel S. Goll<sup>4</sup>, Vanessa Haverd<sup>5,†</sup>, Atul K. Jain<sup>6</sup>, Emilie Joetzjer<sup>7</sup>, Etsushi Kato<sup>8</sup>, Sebastian Lienert<sup>9</sup>, Danica Lombardozzi<sup>10</sup>, Patrick C. McGuire<sup>11</sup>, Joe R. Melton<sup>12</sup>, Julia E. M. S. Nabel<sup>13,17</sup>, Julia Pongratz<sup>13,14</sup>, Stephen Sitch<sup>15</sup>, Anthony P. Walker<sup>16</sup>, and Sönke Zaehle<sup>17</sup>

**Correspondence:** Lina Teckentrup (l.teckentrup@unsw.edu.au)

Received: 10 March 2021 – Discussion started: 15 March 2021

Revised: 8 September 2021 – Accepted: 22 September 2021 – Published: 20 October 2021

**Abstract.** Australia plays an important role in the global terrestrial carbon cycle on inter-annual timescales. While the Australian continent is included in global assessments of the carbon cycle such as the global carbon budget, the performance of dynamic global vegetation models (DGVMs) over Australia has rarely been evaluated. We assessed simulations of net biome production (NBP) and the carbon stored in vegetation between 1901 to 2018 from 13 DGVMs (TRENDY v8 ensemble). We focused our analysis on Australia's short-term (inter-annual) and long-term (decadal to centennial) ter-

restrial carbon dynamics. The TRENDY models simulated differing magnitudes of NBP on inter-annual timescales, and these differences resulted in significant differences in long-term vegetation carbon accumulation (–4.7 to 9.5 Pg C). We compared the TRENDY ensemble to several satellite-derived datasets and showed that the spread in the models’ simulated carbon storage resulted from varying changes in carbon residence time rather than differences in net carbon uptake. Differences in simulated long-term accumulated NBP between models were mostly due to model responses to land-

<footnote n="1"><sup>1</sup>ARC Centre of Excellence for Climate Extremes, Sydney, NSW, Australia</footnote>

<footnote n="2"><sup>2</sup>Climate Change Research Centre, University of New South Wales, Sydney, NSW, Australia</footnote>

<footnote n="3"><sup>3</sup>School of Biological Sciences, University of Bristol, Bristol, BS8 1TQ, UK</footnote>

<footnote n="4"><sup>4</sup>Université Paris Saclay, CEA-CNRS-UVSQ, LSCE/IPSL, Gif sur Yvette, France</footnote>

<footnote n="5"><sup>5</sup>CSIRO Oceans and Atmosphere, G.P.O. Box 1700, Canberra, ACT 2601, Australia</footnote>

<footnote n="6"><sup>6</sup>Department of Atmospheric Sciences, University of Illinois, Urbana, IL 61821, USA</footnote>

<footnote n="7"><sup>7</sup>CNRM, Université de Toulouse, Météo-France, CNRS, Toulouse, France</footnote>

<footnote n="8"><sup>8</sup>Institute of Applied Energy (IAE), Minato-ku, Tokyo 105-0003, Japan</footnote>

<footnote n="9"><sup>9</sup>Climate and Environmental Physics, Physics Institute and Oeschger Centre for Climate Change Research, University of Bern, Bern, Switzerland</footnote>

<footnote n="10"><sup>10</sup>Terrestrial Sciences Section, Climate and Global Dynamics, National Center for Atmospheric Research, Boulder, CO 80305, USA</footnote>

<footnote n="11"><sup>11</sup>Department of Meteorology, Department of Geography and Environmental Science, National Centre for Atmospheric Science, University of Reading, Reading, UK</footnote>

<footnote n="12"><sup>12</sup>Climate Research Division, Environment and Climate Change Canada, Victoria, BC, Canada</footnote>

<footnote n="13"><sup>13</sup>Max Planck Institute for Meteorology, Hamburg, Germany</footnote>

<footnote n="14"><sup>14</sup>Department of Geography, LMU, Munich, Germany</footnote>

<footnote n="15"><sup>15</sup>College of Life and Environmental Sciences, University of Exeter, Exeter EX4 4RJ, UK</footnote>

<footnote n="16"><sup>16</sup>Environmental Sciences Division and Climate Change Science Institute, Oak Ridge National Laboratory, Oak Ridge, Tennessee, USA</footnote>

<footnote n="17"><sup>17</sup>Max Planck Institute for Biogeochemistry, P.O. Box 600164, Hans-Knöll-Str. 10, 07745 Jena, Germany
†deceased, 19 January 2021</footnote>
use change. The DGVMs also simulated different sensitivities to atmospheric carbon dioxide (CO₂) concentration, although notably, the models with nutrient cycles did not simulate the smallest NBP response to CO₂. Our results suggest that a change in the climate forcing did not have a large impact on the carbon cycle on long timescales. However, the inter-annual variability in precipitation drives the year-to-year variability in NBP. We analysed the impact of key modes of climate variability, including the El Niño–Southern Oscillation (ENSO) and the Indian Ocean Dipole (IOD), on NBP. While the DGVMs agreed on sign of the response of NBP to El Niño and La Niña and to positive and negative IOD events, the magnitude of inter-annual variability in NBP differed strongly between models. In addition, we find that differences in the timing of simulated phenology and fire dynamics are associated with differences in simulated or prescribed vegetation cover and process representation. We further find model disagreement in simulated vegetation carbon, phenology, and apparent carbon residence time, indicating that the models have different types and coverage of vegetation across Australia (whether prescribed or emergent). Our study highlights the need to evaluate parameter assumptions and the key processes that drive vegetation dynamics, such as phenology, mortality, and fire, in an Australian context to reduce uncertainty across models.

# 1 Introduction

Decadal variability in the growth rate of atmospheric carbon dioxide (CO₂) is strongly influenced by variability in the uptake and release of carbon by the oceans and the terrestrial biosphere (Ballantyne et al., 2012; Raupach et al., 2008). The inter-annual variability (IAV) in the CO₂ growth rate is dominated by terrestrial processes (e.g. vegetation productivity, respiration, and fire emissions) and their responses to both temperature and precipitation, driven by modes of climate variability (Ahlström et al., 2015; Zhang et al., 2018; Poulter et al., 2014). The El Niño–Southern Oscillation (ENSO) is the dominant mode of global variability (e.g. Ahlström et al., 2015) and contributes around 26% of the IAV in global gross primary production (GPP; Zhang et al., 2019). Globally, ENSO has also been shown to explain more than 40% of satellite-derived net primary production (NPP) variability, mainly driven by the response of Southern Hemisphere ecosystems (Bastos et al., 2013) and in particular semi-arid ecosystems (Zhang et al., 2018).

Around 70% of the land in Australia is classified as either arid or semi-arid (Brown et al., 2008). Oceania has been found to contribute significantly to the uncertainty in global and regional carbon budgets (Bastos et al., 2020), and the important role played by semi-arid ecosystems in explaining the variability in the global carbon cycle was highlighted by the 2011 La Niña event. While on average, ~17% of IAV in net

biome production (NBP) was attributable to Australia over the historical period, during the 2011 La Niña event, Australia contributed around 60% of the global carbon sink (e.g. Poulter et al., 2014). A recent study suggested that the Australian terrestrial carbon sink may be enhanced due to more extreme wet events projected for future decades (Ma et al., 2016). This carbon uptake enhancement has been associated with the asymmetric response of GPP to precipitation (i.e. positive GPP anomalies tend to be larger than negative ones; Haverd et al., 2017) in combination with vegetation expansion linked to rainfall (based on a single dynamic vegetation model; Poulter et al., 2014). At the same time, a series of studies have also identified evidence of rising CO₂, leading to a marked greening of the Australian continent (Donohue et al., 2009, 2013; Ukkola et al., 2016b; Trancoso et al., 2017).

Dynamic global vegetation models (DGVMs) are commonly used to explore large-scale responses of the carbon cycle to climate variability and climate change (Friedlingstein et al., 2019; Le Quéré et al., 2018). However, comparing NBP (the change in carbon stocks including carbon losses due to disturbance) between different DGVMs shows a large model spread with substantial variations in annual global NBP of up to 3 Pg C yr⁻¹ when forced with the same meteorological drivers (Le Quéré et al., 2018). This demonstrates large uncertainties in the representation of the terrestrial carbon cycle and the associated response to climate and land-use change in DGVMs. Roxburgh et al. (2004) concluded from an Australian-based evaluation of NPP that DGVMs were unable to capture the interactions between terrestrial biosphere and atmosphere (models simulated a fivefold variation in annual NPP), citing weaknesses related to the role of water and nutrients in limiting productivity. A number of more recent studies have improved our understanding of the Australian terrestrial carbon cycle, employing flux tower data (e.g. Tarin et al., 2020), satellite-derived carbon fluxes (e.g. Cleverly et al., 2016), and regional biospheric modelling (e.g. Haverd et al., 2013, 2016b).

Nevertheless, while the Australian continent is included in global assessments of the carbon cycle, DGVM performance over Australia has rarely been evaluated. The need for such an evaluation is emphasized by recent studies that have used DGVM simulations to underline the importance of Australia's ecosystems to the inter-annual variability in the global carbon cycle (Poulter et al., 2014; Ahlström et al., 2015). DGVM studies have also previously highlighted vegetation responses to climate as critical to future projected changes in water resources (Ukkola et al., 2016a) and an enhanced future carbon sink (Kelley and Harrison, 2014) across Australia. Thus, evaluating DGVM skill is also equally important for assessing projected changes in Australia's carbon and water cycles.

In this study, we assessed the terrestrial carbon cycle for Australia simulated by 13 DGVMs that are part of the TRENDY v8 ensemble. We examined drivers of both the
long-term (decadal to centennial timescales) and short-term (inter-annual timescales) model responses of the carbon cycle. With these timescales in mind, we expected that model differences result from

- sensitivity to increasing CO₂

- sensitivity to climate variability

- prescribed and simulated land cover and interpretation of land-use change

- assumed and emergent functional representation of Australian vegetation.

We therefore examined the contribution made by each of these differences across the TRENDY v8 ensemble. Our goal was to isolate, to the extent possible, the causes of differences among the TRENDY v8 ensemble as a first step towards prioritizing areas to focus on to resolve the differences among the models.

# 2 Materials and methods

## 2.1 Models and simulations

We used the 13 dynamic vegetation models that are part of the TRENDY v8 model ensemble (Friedlingstein et al., 2019): CABLE-POP, CLASS-CTEM, CLM5.0, ISAM, ISBA-CTRIP, JSBACH, JULES-ES, LPX-Bern, OCN, ORCHIDEE, ORCHIDEE-CNP, SDGVM, and VISIT. In these models, vegetation cover is either dynamically simulated (CABLE-POP, JULES-ES, LPX-Bern, and SDGVM) or prescribed (CLASS-CTEM, CLM5.0, ISAM, ISBA-CTRIP, JSBACH, OCN, ORCHIDEE, ORCHIDEE-CNP, and VISIT). All models used identical forcing inputs and followed a common simulation protocol over the period 1700–2018 (see Fig. C1 and description in the Appendix A for more detailed information). However, each model used a model-specific interpretation of the land-use change forcing.

We remapped all model outputs and satellite datasets (see below) to a common 0.5° grid using first-order conservative regridding (except for the comparison with the data over the North Australian Tropical Transect (NATT)). For this analysis, we defined a year as starting in July and ending in June (except for the C<sub>Veg</sub> analysis, which followed calendar years) to capture the Southern Hemisphere growing season and capture full El Niño and La Niña events, which usually start and end in austral winter. We expressed the change in the variables as the difference to the 1901–1930 average and processed the data with netCDF Operators (NCO; version 4.7.7. http://nco.sf.net, last access: 27 September 2021) and Climate Data Operators (CDO; version 1.9.5. http://mpimet.mpg.de/cdo, last access: 27 September 2021). The data analysis was conducted with Python version 3.

## 2.2 Satellite data

To assess model simulations, we used several satellite-derived datasets. We note that direct continental-scale measurements of the variables below do not exist and that the satellite products also rely on models themselves.

### 2.2.1 Gross primary production and phenology

We compared the simulated GPP phenology cycle to GPP derived from solar-induced chlorophyll fluorescence (SIF) (GOSIF-GPP; Li and Xiao, 2019). To ensure the GOSIF-GPP phenology accurately captured phenology, we also compared these data to a satellite-derived leaf area index (LAI) product (modified Copernicus Service information, 2020; see Fig. C2). GOSIF-GPP is based on the OCO-2-based SIF product (GOSIF) and derived by assuming linear relationships between SIF and GPP estimated from eddy covariances sites. Uncertainties were accounted for by using eight different SIF-GPP relationships (derived universally and biome-specific, with and without intercept at both site and grid cell levels). The 0.05° dataset covers the period 2000 to 2018 at an 8 d temporal resolution that we aggregated to a monthly time step.

### 2.2.2 Carbon stored in vegetation

Estimates of aboveground carbon biomass can be derived from satellite estimates of vegetation optical depth (VOD; related to the aboveground vegetation's water content, density, and biomass). We used aboveground carbon biomass to assess the carbon stored in vegetation in the TRENDY models and note that the datasets are not directly comparable since the satellite product does not account for below-ground carbon. We used global estimates of the annual average aboveground biomass carbon for 1993–2012 with a 0.25° spatial resolution derived from a series of satellite passive microwave instruments (Version 1.0; Liu et al., 2015). Liu et al. (2015) used an empirical approach to convert a harmonized time series of VOD to aboveground biomass carbon.

### 2.2.3 Fire-related CO₂ fluxes

We used two different satellite estimates of fire CO₂ emissions: the Global Fire Emissions Database version 4 (GFED4s) described in van der Werf et al. (2017) and the Copernicus Atmosphere Monitoring Service Global Fire Assimilation System (CMAS-GFAS; Kaiser et al., 2012). GFED4s has a spatial resolution of 0.25° with a monthly time step and provides data from 1997–2020. We note that starting from 2017, the GFED4s data are flagged as a beta version; i.e. starting 2017, fire emissions are derived from active fire detections instead of burned area. CAMS-GFAS covers the years 2003–2020 on a 0.1° horizontal grid. The source for both estimates is the MODIS satellite, but different variables and methods are used to derive the fire CO₂ emissions (for
more detail, see e.g. van der Werf et al., 2017; Kaiser et al., 2012; Li et al., 2019; and Pan et al., 2020).

### 2.2.4 Land-cover fraction

We used the MODIS/Terra Vegetation Continuous Fields Yearly L3 Global 250 m SIN Grid V006 dataset for a comparison of satellite-derived and land cover in DGVMs. This product estimated the global land-cover fraction using surface reflectance, brightness temperature, and the MODIS Global 250 m Land/Water Map (DiMiceli et al., 2017).

## 2.3 Flux sites and North Australian Tropical Transect sites

We analysed the TRENDY models in the context of the NATT gradient, which encompasses a precipitation gradient ranging from ~720 to 1650 mm yr⁻¹. The NATT was established to understand the function of savanna ecosystems, the predominant landscape in northern Australia. We used the four sites Howard Springs (AU-How), Daly Uncleared (AU-DaS), Dry River (AU-Dry), and Sturt Plains (AU-Stp; see Table B1 in Appendix B). We chose the simulations with transient CO₂ concentration, climate, and land-use change. In order to compare the models to the flux sites, we isolated the corresponding grid cells based on the site coordinates and calculated the net ecosystem exchange (NEE) as a balance between GPP and terrestrial ecosystem respiration (TER). We further divided the data into the wet (November–April) and dry (May–October) season.

Although we showed flux tower observations and simulations of NEE, GPP, and TER together, we note that the spatial scales of observed ecosystem fluxes (eddy covariance) and simulations by a DGVM are not directly comparable given that the flux tower footprint is ~1 vs. ~3000 km² for a grid cell in a 0.5° grid. Consequently, DGVMs cannot represent local features such as heterogeneous environmental conditions and land cover or biogeographical feedbacks with atmospheric conditions specific to the site (Piao et al., 2013; Luo, 2007). However, given that the meteorology recorded at the sites and the meteorology used to drive the TRENDY simulations were highly correlated (see Fig. C3 in Appendix C), we assumed the TRENDY models simulated the vegetation in a similar climate as the observed. We also note that due to data loss, some of the observed data contained a high proportion of gap-filled data. This was particularly the case at the Howard Springs site and for the TER flux across sites.

## 2.4 Vegetation classes

To group the analysis by regions, we classified Australia into six vegetation classes (tropics, savanna, warm temperate, cool temperate, Mediterranean, and desert; see Fig. C4) following Haverd et al. (2012). These regions have distinct

climate and biophysical characteristics and are based on the agro-climatic classification by Hutchinson et al. (2005).

## 2.5 Identification of modes of variability

Numerous studies highlighted the impact of different modes of climate variability on Australian weather patterns (e.g. Ummenhofer et al., 2011). The El Niño–Southern Oscillation (ENSO) and the Indian Ocean Dipole (IOD) drive variability in precipitation, which in turn is the main driver of variability in the Australian carbon cycle (e.g. Cleverly et al., 2016; Haverd et al., 2013). We identified years corresponding to the phase of ENSO and IOD events (see Table B2). We used these different climate modes to group our analysis of modelled carbon fluxes.

## 2.6 Apparent carbon residence time in vegetation

To understand model differences that are related to carbon stored in vegetation ($C_{Veg}$), we analysed the carbon residence time in vegetation ($\tau$). We followed Friend et al. (2014) and Pugh et al. (2020), who defined the change in $C_{Veg}$ over time as

$$ \frac{dC_{Veg}}{dt} = NPP - \frac{C_{Veg}}{\tau}. \qquad (1) $$

Consequently, we calculated the change in the carbon residence time using annual time steps according to

$$ \tau = \frac{C_{Veg}}{NPP - \frac{dC_{Veg}}{dt}}. \qquad (2) $$

## 2.7 PFT groups

Each of the TRENDY models has its own vegetation classification. To compare simulated and prescribed vegetation cover between models, we grouped the model-specific plant functional types (PFTs) into 10 PFT groups: evergreen trees (“EVG”); deciduous, summergreen, and raingreen trees (“DCD/SMG/RNG”); shrubs; savanna; $C_3$ grass; $C_4$ grass; $C_3$ agriculture and $C_4$ agriculture (i.e. crops and pasture PFTs); bare ground; and a group comprising the remaining land-cover types, such as city, urban, or lakes. For each model, we assigned the model-specific PFTs to these groups and finally compared the fraction of land covered by each PFT group. Note that not all models account for all the vegetation groups; i.e. only 4 out of the 10 PFT groups included input from all of the models (EVG, DCD/SMG/RNG, $C_3$ grass, and $C_4$ grass).

# 3 Results

## 3.1 Net biome production and carbon stored in vegetation

Summed over Australia, the TRENDY model ensemble simulated large inter-annual variability in NBP (see Fig. 1a)
ranging from −0.9 to 1.4 Pg C yr⁻¹. Figure 1a also shows a considerable model spread: the TRENDY models varied on average by 0.5 Pg C yr⁻¹ with a difference of up to 1.26 Pg C yr⁻¹ between the most extreme models. Years with extremely high or low NBP, such as in 1973 and 2010, were associated with particularly large uncertainty amongst the models. These large differences on inter-annual timescales led to differences in cumulative NBP over 1901–2018, varying between −4.7 (OCN) and 9.5 Pg C (LPX-Bern; see Fig. 1b). Two of the DGVMs simulated a net carbon source over the Australian continent from 1901–2018, while the remaining models simulated a net sink ranging from negligible through to 9.5 Pg C.

Across the 20th Century, the TRENDY models simulated very different amounts of $C_{Veg}$. For example, JSBACH simulated ~2.5 Pg C, while JULES-ES and ORCHIDEE-CNP simulated more than ~16 Pg C (see Fig. 1c). We used an estimate derived from vegetation optical depth (VOD; see Fig. 1c) to put the model simulations into context. A total of 9 of the 13 models predicted a higher amount of $C_{Veg}$ than aboveground biomass estimated by the satellite. Uncertainty in the VOD data was not available, but if we conservatively assume ±30%, at least 7 of the 12 models were outside this range. However, we note that higher $C_{Veg}$ in comparison to aboveground biomass might also result from the fact that $C_{Veg}$ includes above- and belowground biomass. The models also showed very different long-term behaviour. For example, some DGVMs simulated a steady-state $C_{Veg}$ (e.g. JSBACH and LPX-Bern) through 1901–2018, while other models displayed a strong decrease in $C_{Veg}$ until 1970 and a subsequent slight increase (e.g. OCN and ISBA-CTRIP). VISIT simulated a sudden jump in $C_{Veg}$ around 1975. Similarly, $C_{Soil}$ varied strongly between the models, from 11.9 (JSBACH) up to 64.4 Pg C (JULES-ES). Most of the models did not show a change over time, while three models increased in $C_{Soil}$ (CABLE-POP, ISAM, and LPX-Bern). In short, Fig. 1 shows that the 13 TRENDY models simulated different inter-annual variability, contrasting cumulative NBP and inconsistent short- and long-term $C_{Veg}$ and $C_{Soil}$.

## 3.2 Response to atmospheric CO₂ concentration, climate, and land-use change

We used three sets of TRENDY simulations (see Fig. C1) to separate the effect of three drivers – land-use change, atmospheric CO₂ concentration, and climate – on the Australian carbon cycle. Figure 2 shows cumulative NBP over 1901–2018 for the three simulations CO₂, CO₂ + CLIM, and CO₂ + CLIM + LUC. Summed over the years 1901–2018, the TRENDY models showed high variability in the response to the different external drivers. The models agreed on the sign of the response to atmospheric CO₂ and the combined effect of atmospheric CO₂ and varying climate such that, by 2018, all models accumulated NBP on both regional and continen-

tal scales. However, the response of cumulative NBP to increasing atmospheric CO₂ differed strongly among the models, with an increase in NBP ranging from 2.5 up to 10 Pg C over all of Australia. Varying climate in combination with increasing atmospheric CO₂ concentration led to stronger increases in cumulative NBP for some models (e.g. LPX-Bern and VISIT; see Fig. 2; all regions except tropics and Mediterranean) or decreases (e.g. ISAM and ORCHIDEE-CNP). Lastly, land-use change led to the greatest variation among the TRENDY models, with differences in both sign and magnitude (−3.3 to 8.5 Pg C). A total of 8 out of 13 models simulated a decrease in cumulative NBP compared to the CO₂ + CLIM run for all regions. Imposed land use turned positive cumulative NBP in the CO₂ + CLIM simulations into negative cumulative NBP in the CO₂ + CLIM + LUC simulation for OCN and ORCHIDEE-CNP. Only ISAM accumulated more NBP in the CO₂ + CLIM + LUC simulation compared to the CO₂ + CLIM run in most regions.

## 3.3 Response to modes of climate variability

The climate in Australia, especially the precipitation patterns, is strongly influenced by different modes of variability such as ENSO and the IOD. Therefore, the modes of variability can be expected to influence the IAV of the terrestrial carbon cycle over Australia. To further understand the response to the climate forcing, we decomposed the model response to the different modes of variability (see Table B2), specifically to positive and negative IOD events (see “pIOD” and “nIOD” in Fig. 3) and El Niño and La Niña events (see “El Niño” and “La Niña” in Fig. 3). Negative and positive NBP values that are associated with neither ENSO nor the IOD (see “other negative” and “other positive” in Fig. 3) as well as the response of the models to the 10 driest and wettest years (“driest years” and “wettest years”) are also shown in Fig. 3.

In general, the individual models agreed on the sign of NBP for the different climate modes of variability with negative NBP for positive IOD and El Niño events and positive NBP for negative IOD and La Niña events. The two “other” groups (negative and positive) covered a similar value range compared to the modes of variability. For both IOD and ENSO events as well as for the “other negative/positive” categories, the median values across years were similar among the models. Interestingly, the model responses in the driest and wettest years led to the largest range among the models in both median and IAV compared to the other six panels. This highlights the importance of looking beyond climate modes of variability to understand responses of the Australian terrestrial carbon cycle.

## 3.4 Seasonal productivity and phenology

To examine differences in simulated carbon uptake, we evaluated the simulated seasonal cycle. Figure 4 shows that the
<figure><img src="image_2.png" /></figure>

**Figure 1.** Net biome production (NBP), carbon stored in vegetation ($C_{Veg}$), and carbon stored in soil ($C_{Soil}$) for the TRENDY model ensemble for the $CO_2 + CLIM + LUC$ run. Positive values in NBP represent a net carbon sink; negative values are a net carbon source. The solid green line in (a) shows the ensemble mean of total annual NBP, and the shaded green area shows the range across models. Panel (b) shows the cumulative net biome production for each model summed over Australia. Panels (c) and (d) show $C_{Veg}$ and $C_{Soil}$ respectively summed over Australia for each model. The black line in (c) shows a satellite-derived estimate ("VOD"; Liu et al., 2015). Note that we define the first month of the year as July and the last month of the year as June for NBP and follow calendar years for the carbon pools.

TRENDY models varied in the timing and the magnitude of peak productivity for the different vegetation regions. We compared the TRENDY models to the GPP estimate derived from solar-induced chlorophyll fluorescence ("GOSIF-GPP"; Li and Xiao, 2019). We note that because direct measurements of GPP at continental scales do not exist (i.e. GOSIF-GPP is not directly observed GPP), we constrained our comparison to an indicative evaluation of differences. Most models simulated higher productivity for the tropics, savanna, and warm temperate regions compared to GOSIF-GPP. For all regions except the desert, most models were not within the uncertainty band of the GOSIF-GPP satellite estimate (grey band). The timing of peak productivity varied among the TRENDY models and GOSIF-GPP: most models matched peak productivity for the tropics, savanna, and desert. In contrast, for the cool temperate and Mediterranean regions, most models lagged peak productivity estimated by GOSIF-GPP by 1 to 3 months. As the SIF-GPP relationship was derived using a fixed relationship between eddy covariance

and satellite SIF, it is possible that the model data mismatch implies a missing species sensitivity to water stress in the SIF data. However, the consistent timing in peak SIF and satellite LAI (see Fig. C2) tends to imply that the lags in phenology relate to differences in assumed and emergent vegetation cover (see below).

## 3.5 Apparent carbon residence time

Relative to the average over 1901–1930, all models simulated a similar rate in increase in NPP (see Fig. 5a). The apparent carbon residence time ($\tau$), i.e. the balance between growth of plant tissues (leaves, wood, roots) and the turnover of these tissues, varied between 2 and 14 years averaged over Australia. Models that simulated high $C_{Veg}$ (compare Fig. 1) also had high values for $\tau$ (e.g. ORCHIDEE-CNP and JULES-ES) and vice versa (e.g. JSBACH and LPX-Bern). Some models did not simulate changes in $\tau$ over time (e.g. CLM5.0, JSBACH, JULES-ES, LPX-Bern), while other
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Cumulative net biome production (NBP) from 1901–2018. The first group in each panel shows the CO₂ effect on accumulated NBP (“CO₂”; i.e. transient CO₂ forcing, time-invariant climate and land-use mask). The second group (“CO₂ + CLIM”) shows the combined effect of a transient CO₂ and climate forcing (i.e. transient CO₂ and climate forcing, time-invariant land-use mask). The third group (“CO₂ + CLIM + LUC”) shows the combined effect of a transient CO₂, climate forcing, and land-use change (i.e. transient CO₂ and climate forcing and land-use mask). Regions are defined according to Haverd et al. (2012) (see Fig. B4). We define the first month of the year as July, and the last month of the year is defined as June. Note that the y axis limits differ between the panels.</figcaption></figure>

models simulated decreases. In particular, the ISBA-CTRIP and OCN models showed a strong decline – from ~ 6 years (1900) to 4 years (1960) – before levelling off (see Fig. 5c). Even though NPP increased for ISBA-CTRIP and OCN, $C_{Veg}$ significantly declined for these two models as well due to the decrease in $\tau$. The remaining models either balanced increasing NPP and decreasing $\tau$ so that $C_{Veg}$ did not show a strong trend (e.g. CLM5.0, JSBACH, JULES-ES, LPX-Bern) or increased in $C_{Veg}$ because the effect of change in NPP was greater than the decline in $\tau$ (e.g. CLASS-CTEM, ORCHIDEE, VISIT).

## 3.6 Northern Australian Tropical Transect (NATT)

To better understand differences in modelled carbon fluxes, we examined simulations along a rainfall gradient: the Northern Australian Tropical Transect (NATT). Figure 6 shows the probability density functions for NEE, GPP, and TER for the wet (November–April) and dry (May–October) seasons across the NATT. During the dry season, as the sites become drier, the mean position of the GPP and TER distributions shifted closer to zero, and the distribution spread narrows. By contrast, the shape of the distributions varied strongly among the models for all variables: some models (e.g ISAM as well as LPX-Bern for dry season NEE and GPP) simulated high
<figure><img src="image_4.png" /></figure>

**Figure 3.** Net biome production (NBP) during positive IOD ("pIOD"), negative IOD ("nIOD"), El Niño, and La Niña events (see Table B2). All negative NBP values that are associated with neither pIOD nor El Niño are classified as "other negative"; all positive NBP values that are associated with neither nIOD nor La Niña events are classified as "other positive". "Driest years" shows NBP values associated with the 10 driest years, and "wettest years" shows NBP values associated with the 10 wettest years (based on precipitation anomalies averaged over Australia). We define the first month of the year as July, and the last month of the year is defined as June. Note that the y axis limits differ between the panels.

peak probability densities and narrow distributions, while others (e.g. JULES-ES and OCN) displayed a flatter distribution. During the wet season, some models were more productive than observed for all sites (CABLE-POP, JSBACH, JULES-ES, and SDGVM). Overall, there was no clear systematic pattern in the arrangement of the models across the NATT; even when models simulated similar distributions for NEE, they often did so for different reasons (i.e. contrasting trade-offs in carbon uptake and respiration).

## 3.7 Fire-related CO₂ fluxes

Fire is an important component of the carbon cycle in many of Australia's ecosystems. At the monthly timescale, the seven DGVMs that simulated fire dynamics significantly underestimated the fire CO₂ emissions associated with extreme events compared to the GFED4s and CAMS-GFAS satellite products. However, the timing of simulated peak fire events was synchronized with the observations for some models (see Fig. 7a). Fire weather can be associated with ENSO cycles in Australia (e.g. Harris and Lucas, 2019); Fig. 7a however did not indicate a clear connection be-
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Seasonal cycle for satellite-derived gross primary production (GOSIF-GPP; black) and the associated uncertainty (i.e. ±1 SD) as well as the individual TRENDY models averaged over 2000–2017. Regions are defined according to Fig. C4. Note that the y axis limits differ between the panels.</figcaption></figure>

tween ENSO cycles and extreme fire events. The first two peaks in observed fire emissions occur during El Niño events; the two highest peaks in fire CO₂ emission according to CMAS-GFAS instead coincided with the 2011 La Niña event and with the ENSO-neutral year 2012/2013. Figure 7b) shows that all models simulated lower variability compared to the satellite-derived estimates, with similarly wide but less peaked probability density functions. The peak was shifted to the right of the observations for some models (LPX-Bern, SDGVM, and VISIT; see Fig. 7b). On monthly timescales, four TRENDY models capture some features in the variability in the satellite-derived observations with either weakly (ISBA-CTRIP and VISIT: both datasets; JSBACH: GFED4s), moderately (JSBACH: CAMS-GFAS; CLASS-CTEM: CAMS-GFAS), or highly (CLASS-CTEM: GFED4s) significant correlation coefficients. The remaining models do not show a significant relationship to either of the datasets. Aggregated to annual values, the TRENDY models generally underestimated the fire CO₂ emissions and did not

capture the variability in, or timing of, extreme fire years (see Fig. 7c). CLASS-CTEM, JSBACH, and ISBA-CTRIP captured some features of the variability in the satellite-derived observations. CLASS-CTEM is moderately correlated with both datasets ($r > 0.5$), ISBA-CTRIP shows a significant moderate correlation with the GFED4s dataset only, and JSBACH is highly correlated with both datasets ($r > 0.7$; see Table B3). The remaining models are not significantly linked the satellite-derived observations.

## 3.8 Land-cover fraction

Figure 8 shows the vegetation fraction, either prescribed or simulated dynamically, from the TRENDY models averaged over Australia. We grouped the model-specific PFTs into 10 groups (see methods). In general, the average land cover varied quite strongly among the models. The average fraction of Australia covered by natural vegetation differed between 24.9 % (LPX-Bern) and 93.7 % (VISIT). The vegetation composition also displayed large differences, with some
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Change in net primary production (NPP) in comparison to the 1901–1930 average (a), total apparent carbon residence time in vegetation τ (b), change in τ in comparison to the 1901–1930 average (c), and change in carbon stored in vegetation (C<sub>Veg</sub>) in comparison to the 1901–1930 average (d). The solid lines in (b) and (c) show the 5-year moving average; the shaded lines are the original data.</figcaption></figure>

models simulating high tree fractions around 40% (ORCHIDEE and ORCHIDEE-CNP), while other models had low fractions around 2% (e.g. ISAM, and SDGVM). Models accounting for shrubs populated between 3.5% (ISBACTRIP) and 55% (VISIT) of Australia with shrubs. Similarly, land populated by grasses covered a large range between 4% and 50% (VISIT and OCN, respectively), and the fraction of C₃ grasses either exceeded land covered by C₄ grasses or vice versa. Most models showed a relatively small proportion of agricultural land cover except for LPX-Bern. Lastly, models accounting for bare land-cover fraction (i.e. bare ground or desert) defined 0.2%–51% of Australia as bare (VISIT and SDGVM, respectively).

Figure 9 shows the spatial distribution of the fraction of all herbaceous vegetation, the predominant vegetation in Australia, averaged over 1989–2018 for the TRENDY models. The models displayed different patterns in land covered by grass and crops. For example, VISIT had a high grass fraction in almost all regions, with values close to 1, whereas other models showed dense grass cover fractions in the tropics and around the east and west of Australia (e.g. ISAM, ISBA-CTRIP, LPX-Bern, and SDGVM). Some models sim-

ulated a more even distribution of grass cover across Australia with either relatively high (JSBACH) or low fractions (CLASS-CTEM, JULES-ES, and ORCHIDEE-CNP). In contrast, the grass cover fraction for ORCHIDEE generally increased, moving from the coastal areas to the interior of the continent. The tree fraction also varied strongly across the models. Most models had high fractions along the east coast and in the tropics and no significant tree growth in other places, similar to the MODIS satellite-derived tree cover fraction. In contrast, OCN had a low tree cover fraction along the coastline, but this fraction increased towards the interior of Australia, with a higher distribution in the southwest and south-east (see Fig. C5).

## 4 Discussion

On average, Australia's simulated contribution to the global NBP anomaly was 17%, with greater contributions when the global sink (excluding Australia) was small or in years with above-average precipitation in Australia (e.g. in a La Niña years; Poulter et al., 2014). Depending on the model, Aus-
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Probability density functions for observed and modelled net ecosystem exchange (NEE), gross primary production (GPP), and terrestrial ecosystem respiration (TER) for the four flux sites: Howard Springs (AU-How), Daly Uncleared (AU-DaS), Dry River (AU-Dry), and Sturt Plains (AU-Stp). These sites reflect a strong mean annual precipitation (MAP) gradient (provided for each row). We show the monthly distribution of each variable for the wet (November–April) and dry (May–October) seasons. The observed variables originate from eddy covariance data collected by the TERN-OzFlux facility. The probability density function is based on kernel density estimates using Gaussian kernels with bandwidth selection following Scott (1992).</figcaption></figure>

tralia's share of global cumulative NBP could be as high as ~53% (SDGVM). Given Australia's significant contribution to the global terrestrial carbon sink on inter-annual timescales, our goal was to assess the skill of state-of-the-art DGVMs applied to Australia.

Despite being forced with identical meteorology, the 13 DGVMs from the TRENDY v8 ensemble simulated markedly different representations of Australia's terrestrial carbon cycle. The maximum difference in annual NBP among the models reached up to 1.3 Pg C yr⁻¹ (see Fig. 1a). Importantly, uncertainties associated with IAV in simulated NBP accumulated over time, leading to large differences in cumulative NBP among the models by 2018 (see Fig. 1b). The TRENDY models also simulated very different quantities of carbon stocks in vegetation (Fig. 1c). While the satellite estimate of C_Veg lies within the range simulated by the models, individual model estimates of C_Veg varied widely (averaged over 1901–2018 from around 2.7 Pg C for JSBACH to 17.2 Pg C for JULES-ES). Similarly, averaged

over 1901–2018, simulated CSoil differed strongly, from 11.9 (JSBACH) to 64.4 Pg C (JULES-ES).

We identified the key reasons for model differences in NBP related to the processes of land-cover and land-use change, rising atmospheric CO₂ concentration, climate and climate modes of variability, apparent carbon residence time in the vegetation, and fire processes.

## Atmospheric CO₂ concentration

All the models simulated an increase in cumulative NBP in response to increasing CO₂ but with very different magnitudes (+2.46–10.13 Pg C from 1901–2018). This was in line with studies that have identified increasing CO₂ as the main driver of change in the global terrestrial carbon sink and find that model disagreement was largest due to differing sensitivities to CO₂ in cumulative NBP (e.g. Huntzinger et al., 2017; Arora et al., 2020; Walker et al., 2020). The incorporation of nutrient cycles in models and DGVMs in particular tends to result in a reduced sensitivity to CO₂ (e.g. Smith et al., 2014;
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Monthly (a, b) and annual (c) fire CO₂ emissions of seven TRENDY models and two satellite-derived estimates (CMAS-GFAS and GFED4s). Shaded blue areas show La Niña events; shaded red areas show El Niño events. Panels (a) and (c) show the total fire CO₂ emissions; (b) shows the probability density function based on kernel density estimates using Gaussian kernels with bandwidth selection following Scott (1992) of the normalized monthly fire CO₂ emissions. Note that the satellite-derived estimates have a different y axis for the monthly fire CO₂ emissions than for the TRENDY models.</figcaption></figure>

Zaehle, 2013; Thomas et al., 2013; Meiyappan et al., 2015; Meyerholt et al., 2020). In a global analysis, Huntzinger et al. (2017) found that models incorporating a nitrogen cycle tended to have a lower net carbon sink than models without nitrogen constraints. Similarly, two model-based analyses in the context of CO₂ manipulation experiments found lower simulated net primary productivity responses in models that incorporated nitrogen and phosphorus cycles (Medlyn et al., 2016; Fleischer et al., 2019). Interestingly, this finding does not hold for our study: the three models with the highest cumulative NBP all included a nitrogen cycle. A possible explanation may be that Australian soils are considered to be phosphorus- rather than nitrogen-limited (Du et al., 2020; Lambers et al., 2015) due to weathering processes and high phosphorus sorption capacities (Beadle, 1966; Wild, 1958). Results from an ecosystem-scale CO₂ manipulation experiment (Ellsworth et al., 2017) showed that phosphorus availability limited biomass growth in an Australian woodland, although leaf-level photosynthesis was observed to have consistently increased (Jiang et al., 2020; Ellsworth et al., 2017;

Yang et al., 2020). Despite the importance of phosphorus availability in Australia, only two of the TRENDY ensemble models incorporated an interactive phosphorus cycle, and these models (CABLE-POP and ORCHIDEE-CNP) did not show consistent behaviour. It is likely that improved modelling of the phosphorus cycle and in particular the relative “openness” of this cycle and the flexibility of the plant tissue stoichiometry will remain key to accurately simulating the time evolution of Australia’s carbon cycle (Medlyn et al., 2016).

## Climate and modes of climate variability

Variability in Australia's climate can influence the terrestrial carbon sink directly (e.g. fire, droughts or increased water availability; Keenan and Williams, 2018; Ma et al., 2016) or indirectly via interactions with nutrient availability or fire. For example, models that include interactive fire modules might directly (e.g. wind speed driving fire spread, temperature threshold to allow for fire) and indirectly (fuel avail-
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Land-cover fraction in per cent averaged over 1989–2018 over Australia. EVG is the sum of all evergreen tree plant functional types (PFTs); DCD/SMG/RNG is the sum of all deciduous, summergreen, and raingreen tree PFTs; shrubs is the sum of all shrub PFTs; savanna represents the savanna PFT; C₃ grass and C₄ grass include all natural C₃ and C₄ grasses, respectively; C₃ agriculture and C₄ agriculture represent all C₃ and C₄ crops and pasture PFTs, respectively; bare is the sum of bare and desert land-cover types; and other includes land-cover types that are not captured by the above group, such as cities and urban areas as well as lakes. The numbers in brackets represent the land covered by natural vegetation (sum of EVG and RNG, DCD and SMG, shrubs, savanna, and C₃ grasses, and C₄ grasses). Note that not all models account for shrub, savanna, and agriculture PFTs or for explicit bare and desert PFTs.</figcaption></figure>

ability and moisture) depend on climate variables, and altered fire patterns could further affect the simulated terrestrial carbon sink. Overall, we found that the effect of climate on inter-annual NBP was relatively small, with the exception of the Mediterranean region (see Fig. 2). This is in general agreement with Huntzinger et al. (2017), who found that long-term trends in climate were too small to significantly alter the simulated terrestrial carbon sink on long timescales and was a less important factor than CO₂ fertilization, nitrogen deposition (depending on whether the models include an interactive nitrogen cycle), and land-use change. While climate variability might not influence Australian carbon cycle trends on long timescales, inter-annual weather variability, especially the amount and timing of precipitation, is a strong control on inter-annual variability in NBP.

Australia's weather is influenced by different modes of variability, with El Niño and pIOD events tending to result in below-average precipitation in south-eastern Australia (e.g. Ummenhofer et al., 2011). The reduced water availability can be expected to reduce NBP or even turn Australia from a carbon sink to a carbon source (e.g. Ma et al., 2016). In contrast, nIOD and La Niña events are associated with above-average rainfall and consequently an enhanced carbon sink (e.g. Ma et al., 2016; Bastos et al., 2013). The TRENDY

models mostly agreed on sign of the NBP anomaly and simulated negative median values for both pIOD and El Niño events and positive median values for nIOD and La Niña events. The sign of the NBP flux was not unambiguously driven by either ENSO or IOD but was also influenced by other variability (see Fig. 3 “other negative” and “other positive”). The interquartile range for the category “other negative” was similar to those for the pIOD and El Niño categories and displayed more negative outliers than the other modes of variability. This implies that years with extremely low NBP were not necessarily associated with pIOD or El Niño but instead driven by other modes of climate variability and/or periods that did not reach the threshold to be defined as an ENSO or IOD event. Especially for the nIOD, La Niña, and other positive categories, the models varied strongly in the interquartile range, indicating that the individual models simulated different responses to the modes of variability. Studies showed that other modes of variability, such as the Southern Annual Mode (SAM), can influence the Australian weather patterns and consequently the terrestrial vegetation (e.g. Cleverly et al., 2016). We did not include an analysis of SAM because these events are short-lived (1–2 weeks), but future work may examine the resulting impacts. We note that selecting years based on ENSO and IOD events does not
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Land-cover fraction of non-tree vegetation averaged over 2001–2018. Non-tree vegetation here shows the sum of all natural C₃ and C₄ grasses as well as C₃ and C₄ pasture and agriculture and savanna PFTs in the TRENDY models. The MODIS panel shows the non-tree fraction variable from the MODIS/Terra Vegetation Continuous Fields dataset. White areas are missing values.</figcaption></figure>

completely isolate the effect of ENSO or IOD events on the terrestrial carbon cycle as these modes of variability are not independent and could have legacy impacts that last more than 1 year. Further, years where both an IOD and an ENSO event occur were accounted for in both the IOD and ENSO categories and therefore double-counted. An additional limitation of our approach was that we compiled the occurrence of ENSO and IOD events from different sources based on the sea surface temperature. Observations of sea surface temperature prior to 1960 however are associated with large uncertainty given that they were mostly based on ship data (e.g. Deser et al., 2010). In addition, the Southern Hemisphere tends to be less well observed, and spurious trends in analyses do occur (e.g. Hines et al., 2000), leading to an overall increased uncertainty in the sampling of data based on modes of variability.

Overall, dry and wet extremes in precipitation led to the largest differences among the models (see “driest years” and “wettest years” in Fig. 3), suggesting that the models’ sensitivities to wet and dry extremes vary strongly. About half of the TRENDY ensemble had a stronger response to wet than

to dry extremes. This is in line with studies that have noted an increased carbon uptake resulting from the asymmetry in the inter-annual distribution of rainfall as well as the asymmetry in the response of GPP to precipitation (e.g. Haverd et al., 2017). The ecosystem in Australia largely consists of vegetation that is adapted to drought (e.g. Cleverly et al., 2016; Li et al., 2018; De Kauwe et al., 2020). About half of the TRENDY ensemble displayed greater divergence during dry extremes, implying that models did not accurately capture process responses as water becomes limiting. This result is consistent with a model intercomparison at an Australian woodland site, which showed that disagreement was greatest in low-rainfall years (Medlyn et al., 2016). The TRENDY models simulated not only different carbon uptake sensitivities to precipitation but also different seasonal phenology (see Fig. 4), suggesting that differing model sensitivities to rainfall may result as much from the underlying simulated vegetation as the different mechanisms at play in wet and dry extremes. We further looked at the TRENDY output in the context of the NATT, defined by a strong rainfall gradient (~720 to 1650 mm yr⁻¹). We found that for dry-season
GPP, the mean position of the distributions shifts closer to zero, and the distribution narrows as the sites become more arid, implying that vegetation is only productive following rain (Whitley et al., 2016). By contrast, there was no clear systematic pattern in the responsiveness of simulated productivity during the wet and dry seasons.

## Apparent carbon residence time in vegetation

The capacity of the terrestrial vegetation to store carbon depends on the magnitude of the input carbon flux and the carbon residence time in plant tissues (τ; Luo et al., 2003). All models simulated increases in NPP at similar rates compared to their initial NPP (1901–1930 average; see Fig. 5). As the other models simulated a similar change in the rate of carbon uptake throughout the period, divergence in the change in $C_{Veg}$ instead resulted from differences in the change in τ. τ depends on the turnover of plant tissues, the simulated mortality rates, and mortality induced by competition. Each of these processes can be affected indirectly through shifts in vegetation composition (Friend et al., 2014). Previous studies have found that τ was a key source of model uncertainty shared amongst DGVMs (Friend et al., 2014; Pugh et al., 2020). In an analysis of CMIP5 models, Carvalhais et al. (2014) found that fewer than a quarter of the models were within the range of observed τ in tropical Australia and parts of the north-eastern savannas. Similarly, we found that τ in the TRENDY models was associated with large uncertainty: between 2 and up to 14 years over Australia. Models that simulated high vegetation carbon storage also had comparably long τ and vice versa. These differences in modelled τ imply that the TRENDY models assume different parametrizations associated with their PFTs, mortality, and/or simulated different vegetation compositions, and this is an area requiring future evaluation.

## Fire-related CO₂ fluxes

The seven TRENDY models that included fire outputs consistently underestimated satellite-derived peak monthly emissions and did not accurately capture the timing of peak fire CO₂ emissions. Importantly, we found that DGVMs did not always accurately capture fire dynamics in Australia, even when aggregated annually (which should reduce timing biases). While land-use change has been identified to be a significant driver of global burned area (Teckentrup et al., 2019), we found that the effect of land-use change on fire CO₂ emissions in Australia was small (see Fig. C7). Models either increase or decrease in fire CO₂ emissions due to land-use change, but these changes were mostly smaller than 1%. All models simulated fire with varying degrees of complexity and derive fire CO₂ emissions based on model-specific burned area using emission factors based on Andreae and Merlet (2001) and Akagi et al. (2011). Consequently, model

disagreements must originate from differences in simulated burned area.

Biases in burned area have been found to be associated with underlying vegetation composition (which determines fuel load; Teckentrup et al., 2019) and how vegetation cover changes with rising CO₂ and land-use change (through the conversion of natural vegetation to agriculture; Teckentrup et al., 2019). Model evaluations suggest that fire models broadly capture first-order patterns of emissions under present-day conditions (e.g. Li et al., 2019; Hantson et al., 2020), but simulations of seasonality and inter-annual variability are associated with high uncertainty due to differences in process representation. We did not find evidence that models linked extreme fire events with ENSO as suggested by Harris and Lucas (2019). However, the limited availability of “observed” fire CO₂ emissions makes it difficult to draw firm conclusions. In addition, while Harris and Lucas (2019) found that El Niño-like conditions can be associated with extreme fire weather, they also emphasized that this depends on the regions and is most pronounced in eastern Australia, whereas our analysis focused on all of Australia. Lastly, global fire emission estimates based upon satellite observations are still associated with substantial uncertainty, as reflected by the considerable differences in spatial and temporal patterns between different data products (e.g. Li et al., 2019). To account for satellite-derived uncertainty, we included two satellite products, the GFED4s (van der Werf et al., 2017) and the CAMS-GFAS (Kaiser et al., 2012) datasets, that are derived based on different variables (burned area until 2016 and subsequently active fire detections for GFED4s, fire radiate power for CAMS-GFAS). However, since both products rely on the MODIS sensor, they are not independent and consequently only capture the uncertainty associated with satellite-derived fire CO₂ emission estimates to a limited extent.

## Land-cover and land-use change

We found that the models either prescribed or simulated very different fractions of woody and herbaceous cover (see Figs. 9 and C5) and that these discrepancies likely explained a large proportion of the model divergence in simulated $C_{Veg}$ and τ (see below). Notably, the two models with the lowest $C_{Veg}$ had a comparably low tree cover (JSBACH and LPX-Bern; see Fig. 8) and a high proportion of bare land. In contrast, the two models with the highest $C_{Veg}$ (JULES-ES and ORCHIDEE-CNP) simulated a high tree cover fraction. However, Fig. 8 also implies that seemingly similar PFTs were parametrized in different ways in the individual models. For example, CLM5.0 and JSBACH simulated a similar fraction of land covered by natural vegetation, with a comparable vegetation composition comprising trees, shrubs, natural grass, agriculture, and bare land. However, CLM5.0 stored around 4 times more carbon in vegetation compared to JSBACH.
Land-use and land-cover change is one of the major
drivers of the terrestrial carbon sink on global and regional
scales (e.g. Pugh et al., 2019; Huntzinger et al., 2017). Over
Australia, the largest divergence in cumulative NBP was as-
sociated with the implementation of land-use change (LUC;
see Fig. 2). Overall, the TRENDY models simulated different
magnitudes of cumulative NBP in response to LUC, ranging
from -15.5 (OCN) up to 2.9 PgC (ISBA-CTRIP). For two
models, accounting for LUC turned almost all regions from
a carbon sink to a carbon source (OCN, ORCHIDEE-CNP).

Whilst DGVMs were developed to understand the inter-
actions between natural vegetation and the atmosphere, the
insight that land-use is a major driver of the terrestrial car-
bon sink has driven efforts to incorporate land-use change
into DGVMs (Marquer et al., 2018; Pongratz et al., 2018).
Differences in the models' response to LUC result from pro-
cess inclusions, parametrizations, and the model-specific in-
terpretation of the land-use forcing (e.g. change in vegetation
composition as well as wood and crop harvest). Although
LUC information was taken from the Land-Use Harmoniza-
tion 2 (LUH2) guidelines, the modellers had to decide how
to implement it. For example, the LUH2 suggests that all
natural vegetation should be cleared for the conversion to
managed pasture, while the conversion to rangeland only re-
quires the transformation of forest to (managed) grassland.
Other natural vegetation, such as grass- or shrubland, was
used for animal browsing without any transformation of the
land-cover type (Ma et al., 2020). Furthermore, the fraction
of croplands reconstructed by LUH2 increased markedly in
south-western Australia and south-eastern Australia. This in-
crease in crop fraction was associated with increased harvests
and consequently expected to reduce NBP. Some models ex-
hibited this feature to varying degrees (JULES-ES, OCN,
and ORCHIDEE-CNP), while other models barely increased
the crop fraction and impact on NBP (CABLE-POP, OR-
CHIDEE, and VISIT; see Fig. C6). These different represen-
tations of LUC led to emergent differences (e.g. in physiol-
ogy as the models transition between PFTs and agricultural
management) when harvesting was imposed on model simu-
lations.

5 Future directions

Our analysis highlights that considerable work is still re-
quired to improve our capacity to simulate the dynamics of
Australian vegetation and the resulting terrestrial carbon cy-
cle. This finding is consistent with previous model evalua-
tions of NPP (Roxburgh et al., 2004), carbon flux responses
to rainfall (Whitley et al., 2016, 2017), and responses to ele-
vated CO2 (Medlyn et al., 2016).

When DGVMs are applied to Australia, the assumed plant
traits are not adjusted to reflect differences in observed phys-
iology (e.g. Bloomfield et al., 2019; Cernusak et al., 2011),
investment in plant tissues (e.g. Togashi et al., 2015), how

moisture limits plant function (e.g. De Kauwe et al., 2020; Li et al., 2018), or responses of photosynthesis and growth to temperature (Drake et al., 2015). For example, models typically assume much lower values for the key model parameter Vcmax (the maximum rate of carboxylation by the enzyme rubisco; Rogers, 2013) than is commonly observed in Australia (e.g. Bloomfield et al., 2019; Cernusak et al., 2011), implying that DGVMs miss key carbon and water trade-offs. Whitley et al. (2017) previously stressed the importance of improving process representation of phenology, root-water uptake, rooting depth, and fire to improve the representation of savanna ecosystem dynamics in DGVMs. These savanna ecosystems cover around 20% of the Australian landscape, but as this study implies (see Fig. 6), progress has been limited, although recent model developments incorporating competing rooting strategies and dynamic root growth (Sakschewski et al., 2020) may help this process. Australia has a high fraction of endemic vegetation (> 90%; Chapman, 2006), structurally distinct vegetation including open-canopy forests and woodlands (predominately Eucalyptus) and hummock grasslands, and a dominance of sclerophyll leaves (Box, 1996). Thus, there remains a significant opportunity in future work to directly test the current hypothesis implicitly embedded in DGVMs, i.e. that Australia's vegetation is not distinct from other continents. The forthcoming AusTraits (https://ardc.edu.au/project/austraits-accurated-plant-trait-database-for-the-australian-flora/, last access: 8 March 2021) database presents one possibility to examine how the use of Australia-specific plant traits may change DGVM simulations.

One core area in which model development should be
focused relates to the processes that govern the vegetation
cover. In these models vegetation cover is either dynami-
cally simulated (CABLE-POP, JULES-ES, LPX-Bern, and
SDGVM) or prescribed (CLASS-CTEM, CLM5.0, ISAM,
ISBA-CTRIP, JSBACH, OCN, ORCHIDEE, ORCHIDEE-
CNP, and VISIT). Our study has highlighted major differ-
ences across models under historic forcing (see Figs. 8 and
9), implying that future simulations need to be interpreted
cautiously. Haxeltine et al. (1996) demonstrated a realistic
simulation of Australia's potential vegetation cover with the
BIOME2 model, and yet, 25 years later, it is not clear that
DGVMs have made significant progress. One explanation for
this apparent discrepancy may simply relate to the focus of
the modellers. In the Haxeltine et al. (1996) study, the model
explicitly aimed to reproduce Australia's potential vegetation
cover, whereas the current TRENDY models are instead fo-
cused on global applications.

Nevertheless, our results imply that a one-size-fits-
all approach does not appear to work for Australia.
Australia's National Vegetation Information System
(https://www/environment.gov.au/land/native-vegetation/
national-vegetation-information-system, last access:
8 March 2021) generates maps of both vegetation dis-
tribution and structure, which could be used to evaluate
and develop DGVM predictions of potential vegetation. However, a direct comparison to observed vegetation distributions still requires model-specific rules to transform simulated foliar projected cover and height to vegetation classes. Our results also showed that imposing land-use change adds further uncertainty to the vegetation distribution. Despite the importance of improving simulated vegetation type, as Fisher et al. (2015) highlighted, it is not necessarily a straightforward task, but overreliance on bioclimatic boundaries is not preferable moving forwards. Emergent differences in cover type may also originate from parameter assumptions, which determine carbon uptake and loss, additional hypotheses that govern demography and responses to fire. Haverd et al. (2016a) demonstrated promise in accurately predicting tree-grass partitioning in Australian savannas using optimality theory to determine daily carbon allocation. However, this approach has yet to be incorporated into any of the Australian land models (e.g. CABLE-POP), so a fuller evaluation is still required, including the resulting interactions with a nutrient cycle. Kelley et al. (2014) demonstrated significant progress in both the simulation of fires in Australia and the resulting grass-tree cover fractions by incorporating PFT-specific fire (adaptive bark thickness) resistance and post-fire epicormic resprouting into a DGVM. However, these features have yet to be incorporated into any current DGVMs when applied to Australia.

In conclusion, there is a major opportunity in evaluating DGVMs using Australia as a laboratory. The strong inter-annual variability in precipitation offers an opportunity to evaluate the simulated response of the Australian ecosystems to both extreme wet events and droughts. Similarly the repeated frequency of heat extremes (van der Horst et al., 2019; van Gorsel et al., 2016; Mitchell et al., 2014) offers important tests of underlying physiology assumptions. It is likely that insights gained from model evaluation in Australia's environment of longer-duration and high-intensity climate extremes, including droughts, heatwaves, and floods, will help understand how Northern Hemisphere systems will respond to extremes and in particular events that may intensify in the future (IPCC, 2012). Finally, our study focused on Australia's carbon cycle, but the carbon and water cycle are intimately linked, and divergence between model simulations of the carbon cycle implies further work may also be needed to evaluate DGVM simulation of the hydrological cycle.

# Appendix A: Forcing data and simulation protocol

## A1 Climate forcing

The modelling groups chose either the 0.5° Climatic Research Unit (CRU) monthly historical forcing over 1901–2018 (Harris, 2019a) or the 0.5° CRU-JRA55 6-hourly historical forcing over 1901–2018 (Harris, 2019b) regridded to

the CRU 0.5° grid. The variables temperature, downward solar radiation flux, specific humidity, and precipitation in JRA-55 are aligned to temperature, cloud fraction, vapour pressure, and precipitation in CRU TS (v4.03), respectively. Atmospheric pressure, downward long-wave radiation flux, and the meridional and zonal components of wind speed are not modified. JRA-55 is used for the years from 1958 to 2018. For years between 1901 and 1957, random (but fixed) years from JRA-55 for 1958–1967 are used to fill.

## A2 CO₂ concentration

The atmospheric CO₂ concentration is derived from ice core CO₂ data merged with NOAA data from 1958 onwards. The forcing covers the years 1700–2018 incremented annually (Le Quéré et al., 2018). The data from March 1958 are monthly averages from the Mauna Loa (MLO) and the South Pole Observatory (SPO) provided by NOAA's Earth System Research Laboratory (http://www.esrl.noaa.gov/gmd/ccgg/trends/, last access: 8 March 2021). When no SPO data are available (including prior to 1975), SPO is constructed from the 1976–2017 average MLO–SPO trend and average monthly departure. Data for 2016–2018 are preliminary values. Data prior to March 1958 are estimated with a cubic spline fit to ice core data from Joos and Spahni (2008).

## A3 Land use

The land-use datasets are based on updated data from HYDE for the years 1960–2019 as well as the latest wood harvest data provided by the Food and Agriculture Organization (FAO). The land-use states and transitions are identical to the LUH2 v2h dataset (Hurtt et al., 2020) for the years from 1700 to 1949 (i.e. consistent with the input for CMIP6). Starting in 1950, the land-use forcing is based on new inputs from HYDE and new FAO data for the national wood harvest demands. This leads to differences in comparison to LUH2 v2h, primarily in Brazil.

In order to convert natural vegetation to managed pasture, the LUH2 guidelines suggested all natural vegetation is cleared. The conversion of natural vegetation to range-land only requires the clearance of forests. However, each modelling group developed model-specific land-use forcings leading to inter-model differences. We show the change in agricultural land cover in Fig. C6.

## A4 Nitrogen deposition and fertilization

For the years from 1850 to 2014, the TRENDY models were forced with the historical nitrogen deposition database and then transitioned to the future RCP8.5 nitrogen deposition databases for the years from 2015 to 2018. Details are available at https://esgf-node.llnl.gov/search/input4mips/. Nitrogen fertilizer input datasets are available via the N₂O Model Intercomparison Project (“NMIP”; Tian et al., 2018). NMIP assumes that the nitrogen input data remain unchanged in the
years 2015, 2016, 2017, and 2018. Manure is not included for the TRENDY simulations.

## A5 Simulation protocol

The TRENDY models either simulated or prescribed vegetation cover. All models prescribed land-use change according to the LUH2 guidelines. All models used the atmospheric CO₂ concentration from the year 1700 (276.6 ppm) and recycled the climate mean and variability from the years 1901–1920 for the spin-up. Land-use change was constant during the spin-up and the crops, and pasture distribution was set constant to the 1700 set-up. Depending on the simulation, different drivers became transient starting in 1701.

We use three simulations of the TRENDY experiment (see Fig. C1; Friedlingstein et al., 2019): the run with transient atmospheric CO₂ concentration and time-invariant climate forcing and land-use (“CO₂”); the run with transient atmospheric CO₂ concentration and climate forcing and time-invariant land-use (“CO₂ + CLIM”); and finally the run with transient atmospheric CO₂ concentration, climate forcing, and land-use change (“CO₂ + CLIM + LUC”).

### A5.1 Spin-up

For the spin-up, the atmospheric CO₂ concentration and land cover were set to the pre-industrial values of the year 1700. During the spin-up, the years from 1901 to 1920 from the climate forcing were recycled.

### A5.2 CTRL

After the spin-up, the atmospheric CO₂ concentration remained time-invariant for the years 1701–2018. The simulation continued to recycle the 1901–1920 climate forcing, and the land-use was prescribed to the distribution of the year 1700.

### A5.3 CO₂

After the spin-up, the CO₂ run used a transient atmospheric CO₂ forcing. The simulation continued to recycle the 1901–1920 climate forcing, and the land-use was prescribed to the distribution of the year 1700 (as in the “CTRL” run).

### A5.4 CO₂ + CLIM

After the spin-up, the CO₂ + CLIM run used a transient atmospheric CO₂ forcing. The simulation continued to recycle the 1901–1920 climate forcing until the year 1900, and the land-use was prescribed to the distribution of the year 1700. Starting from 1901, the climate forcing became transient.

### A5.5 CO₂ + CLIM + LUC

After the spin-up, the CO₂ + CLIM run used a transient atmospheric CO₂ and land-use forcing, while the simulation continued to recycle the 1901–1920 climate forcing until the year 1900. Starting from 1901, the climate forcing became transient as well.

The TRENDY protocol required all participating models to be in equilibrium for the CTRL simulation after the spin-up. Further, all models had to simulate the net annual land flux as a carbon sink over the 1990s and/or 2000s for the CO₂ + CLIM + LUC run as constrained by global atmospheric and oceanic observations (Keeling and Manning, 2014). Lastly, the global net annual land-use flux (ELUC) had to be a carbon source over the 1990s (based on the CO₂ + CLIM and CO₂ + CLIM + LUC simulation).
Appendix B: Tables

<table><caption>**Table B1.** Information for the four flux sites located in the Northern Australian Tropical Transect (NATT). This includes the mean annual precipitation (MAP), mean annual temperature (MAT), and the vegetation classification following the International Geosphere-Biosphere Programme (IGBP). The observed variables originate from eddy covariance data collected by the TERN-OzFlux facility. Further details of the site vegetation are available in Whitley et al. (2016).</caption><thead><tr><th></th><th>Howard Springs<br>(AU-How)</th><th>Daly Uncleared<br>(AU-DaS)</th><th>Dry River<br>(AU-Dry)</th><th>Sturt Plains<br>(AU-Stp)</th></tr></thead><tbody><tr><td>Years</td><td>2002–2018<br>(2002–2019)</td><td>2008–2019<br>(2008–2018)</td><td>2009–2018<br>(2008–2019)</td><td>2009–2018<br>(2008–2018)</td></tr><tr><td>Latitude</td><td>−12.4952</td><td>−14.1592</td><td>−15.2588</td><td>−17.1507</td></tr><tr><td>Longitude</td><td>131.15005</td><td>131.3881</td><td>132.3706</td><td>133.3502</td></tr><tr><td>MAP [mm yr<sup>−1</sup>]</td><td>1662</td><td>1324</td><td>925</td><td>719</td></tr><tr><td>MAT [°C]</td><td>26.8</td><td>26.8</td><td>27.0</td><td>26.3</td></tr><tr><td>Vegetation</td><td>Woody Savannas:<br>herbaceous and other<br>understorey systems,<br>30 %–60 % forest canopy</td><td>Savannas: herbaceous<br>and other understorey<br>systems, 30 %–60 %<br>forest canopy</td><td>Savannas: herbaceous<br>and other understorey<br>systems, 30 %–60 %<br>forest canopy</td><td>Grasslands: herbaceous<br>types of cover; tree and<br>shrub cover is less than<br>10 %</td></tr></tbody></table>

<table><caption>**Table B2.** Years from 1901–2017 identified as El Niño, La Niña, positive Indian Ocean Dipole (“pIOD”), negative Indian Ocean Dipole event (“nIOD”), and remaining years that are both ENSO- and IOD-neutral based on the analysis by the NOAA (for the ENSO events) and Bureau of Meteorology (IOD events). Note that we define a year as starting in July and ending in June so that for example the 1991–1992 El Niño appears in this table as the 1991 El Niño. We follow Ummenhofer et al. (2011), the Bureau of Meteorology Australia (Bureau of Meteorology, Commonwealth of Australia, 2016), and the National Weather Service Climate Prediction Center (Climate Prediction Center Internet Team, 2021) for the identification of ENSO and IOD events.</caption><thead><tr><th colspan="2">El Niño</th><th colspan="2">La Niña</th><th colspan="2">pIOD</th><th colspan="2">nIOD</th><th colspan="2">ENSO and IOD<br>neutral</th></tr></thead><tbody><tr><td>1902</td><td>1972</td><td>1903</td><td>1964</td><td>1902</td><td>1963</td><td>1906</td><td>1974</td><td>1901</td><td>1951</td></tr><tr><td>1905</td><td>1976</td><td>1906</td><td>1970</td><td>1913</td><td>1972</td><td>1909</td><td>1975</td><td>1904</td><td>1952</td></tr><tr><td>1911</td><td>1977</td><td>1909</td><td>1973</td><td>1919</td><td>1982</td><td>1915</td><td>1980</td><td>1907</td><td>1953</td></tr><tr><td>1914</td><td>1982</td><td>1910</td><td>1974</td><td>1923</td><td>1983</td><td>1916</td><td>1981</td><td>1908</td><td>1959</td></tr><tr><td>1923</td><td>1986</td><td>1916</td><td>1975</td><td>1926</td><td>1995</td><td>1917</td><td>1985</td><td>1912</td><td>1962</td></tr><tr><td>1925</td><td>1987</td><td>1917</td><td>1988</td><td>1935</td><td>1997</td><td>1930</td><td>1989</td><td>1920</td><td>1966</td></tr><tr><td>1930</td><td>1991</td><td>1922</td><td>1989</td><td>1944</td><td>2002</td><td>1933</td><td>1992</td><td>1921</td><td>1967</td></tr><tr><td>1940</td><td>1994</td><td>1924</td><td>1993</td><td>1945</td><td>2006</td><td>1942</td><td>1996</td><td>1927</td><td>1971</td></tr><tr><td>1941</td><td>1997</td><td>1928</td><td>1998</td><td>1946</td><td>2012</td><td>1958</td><td>1998</td><td>1929</td><td>1978</td></tr><tr><td>1957</td><td>2002</td><td>1933</td><td>1999</td><td>1957</td><td>2015</td><td>1960</td><td>2010</td><td>1931</td><td>1979</td></tr><tr><td>1963</td><td>2004</td><td>1938</td><td>2000</td><td>1961</td><td></td><td>1964</td><td>2014</td><td>1932</td><td>1984</td></tr><tr><td>1965</td><td>2006</td><td>1942</td><td>2007</td><td></td><td></td><td></td><td></td><td>1934</td><td>1990</td></tr><tr><td>1968</td><td>2009</td><td>1949</td><td>2008</td><td></td><td></td><td></td><td></td><td>1936</td><td>2001</td></tr><tr><td>1969</td><td>2015</td><td>1950</td><td>2010</td><td></td><td></td><td></td><td></td><td>1937</td><td>2003</td></tr><tr><td></td><td>1954</td><td>2011</td><td></td><td></td><td></td><td>1939</td><td></td><td>2005</td><td></td></tr><tr><td></td><td>1955</td><td>2012</td><td></td><td></td><td></td><td>1943</td><td></td><td>2013</td><td></td></tr><tr><td></td><td>1956</td><td></td><td></td><td></td><td></td><td>1947</td><td></td><td>2016</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td>1948</td><td></td><td></td><td></td></tr></tbody></table>
<table><caption>**Table B3.** Fire CO₂ emissions for Australia averaged over 2003–2018 and the Pearson correlation coefficients between the TRENDY models and the respective observation data. Bold numbers indicate a significant statistical relationship (p value < 0.05).</caption><thead><tr><th rowspan="2">Model</th><th rowspan="2">Fire CO<sub>2</sub> emissions<br>(Pg C yr<sup>-1</sup>)</th><th colspan="2">R (CAMS-GFAS,<br>model)</th><th colspan="2">R (GFED4s,<br>model)</th></tr><tr><th>Monthly</th><th>Annual</th><th>Monthly</th><th>Annual</th></tr></thead><tbody><tr><td>CLASS-CTEM</td><td>0.010 ± 0.006</td><td><strong>0.68</strong></td><td><strong>0.55</strong></td><td><strong>0.78</strong></td><td><strong>0.61</strong></td></tr><tr><td>CLM5.0</td><td>0.014 ± 0.013</td><td>0.12</td><td>0.21</td><td>0.08</td><td>0.27</td></tr><tr><td>ISBA-CTRIP</td><td>0.028 ± 0.021</td><td><strong>0.35</strong></td><td>0.50</td><td><strong>0.41</strong></td><td><strong>0.53</strong></td></tr><tr><td>JSBACH</td><td>0.012 ± 0.007</td><td><strong>0.50</strong></td><td><strong>0.84</strong></td><td><strong>0.43</strong></td><td><strong>0.84</strong></td></tr><tr><td>LPX-Bern</td><td>0.006 ± 0.002</td><td>0.05</td><td>0.41</td><td>0.04</td><td>0.39</td></tr><tr><td>SDGVM</td><td>0.048 ± 0.011</td><td>0.10</td><td>0.14</td><td>0.09</td><td>0.20</td></tr><tr><td>VISIT</td><td>0.026 ± 0.010</td><td><strong>0.21</strong></td><td>0.17</td><td><strong>0.26</strong></td><td>0.18</td></tr><tr><td>CAMS-GFAS</td><td>0.043 ± 0.049</td><td></td><td></td><td></td><td></td></tr><tr><td>GFED4s</td><td>0.030 ± 0.029</td><td></td><td></td><td></td><td></td></tr></tbody></table>

## Appendix C: Figures

<figure><img src="image_11.png" /><figcaption>**Figure C1.** Simulations conducted by the TRENDY models (Friedlingstein et al., 2019). Italic forcing datasets are not used by all models.</figcaption></figure>
<figure><img src="image_12.png" /></figure>

**Figure C2.** Seasonal cycle for satellite-derived gross primary production (GOSIF-GPP; black) and satellite-derived leaf area index data (LAI; Copernicus Global Land Service; green). Shaded areas indicate the according uncertainty (i.e. ± 1 SD.).
<figure><img src="image_13.png" /><figcaption>**Figure C3.** Monthly temperature, precipitation, and incoming short-wave radiation data that are observed (blue) and the corresponding data from the reanalysis product CRU-JRA (orange) for the four flux sites located in the Northern Australian Tropical Transect (NATT): Howard Springs (AU-How), Daly Uncleared (AU-DaS), Dry River (AU-Dry), and Sturt Plains (AU-Stp). The figure titles display the Pearson correlation coefficient between observation and reanalysis ($\rho_{mon}$). The bottom panels show observed monthly NEE (black) and simulated NEE by the individual TRENDY models.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure C4.** Different vegetation classes in Australia according to Haverd et al. (2012).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure C5.** Land-cover fraction of tree vegetation averaged over 2001–2018. Tree vegetation here shows the sum of tree, forest, and shrub PFTs in the TRENDY models. The MODIS panel shows the tree fraction variable from the MODIS/Terra Vegetation Continuous Fields dataset. White areas are missing values.</figcaption></figure>
<figure><img src="image_16.png" /></figure>

**Figure C6.** Cumulative NBP for all TRENDY models over 1901–2018. Land-use change based on Hurtt et al. (2020): increase in C₃ and C₄ crops, rangeland, managed pasture, and urban areas in 2018 compared to 1901.

<figure><img src="image_17.png" /></figure>

**Figure C7.** Difference between CO₂ + CLIM + LUC and CO₂ + CLIM simulations for monthly (a) and annual (b) fire CO₂ emissions of six TRENDY models.
**Code and data availability.** All eddy covariance data are available from http://www.ozflux.org.au/ (OzFlux, 2017). The Global Above-ground Biomass Carbon (version 1.0) dataset (Liu et al., 2015) is freely available from http://wald.anu.edu.au/global-biomass/, and the GOSIF-GPP product (Li and Xiao, 2019) can be obtained from https://globalecology.unh.edu/data/GOSIF-GPP.html. Fire CO₂ emissions were provided by the Copernicus Atmosphere Monitoring Service Global Fire Assimilation System (CAMS GFAS; https://apps.ecmwf.int/datasets/data/cams-gfas/, ECMWF, 2021), and the Global Fire Emissions Database version 4 (GFED4s) described in van der Werf et al. (2017) is available from https://globalfiredata.org/pages/data/. The MODIS/Terra Vegetation Continuous Fields dataset was provided by NASA's Land Processes Distributed Active Archive Center (https://lpdaac.usgs.gov/products/mod44bv006/, DiMiceli et al., 2017). The TRENDY v8 model output is available upon request (https://sites.exeter.ac.uk/trendy, last access: 27 September 2021). All analysis scripts are accessible on https://doi.org/10.5281/zenodo.5570974 (Teckentrup, 2021).

**Author contributions.** LT, MGDK, and AJP designed the experimental analysis. LT wrote the code and analysed the TRENDY models. DSG, VH, AKJ, EJ, EK, SL, DL, PCM, JRM, JEMSN, and SZ conducted the model runs to produce the TRENDY data. LT wrote the first draft with contributions from MGDK and AJP. All authors contributed to the final manuscript. Correspondence and requests for materials should be addressed to LT (l.teckentrup@unsw.edu.au).

**Competing interests.** Some authors are members of the editorial board of *Biogeosciences*. The peer-review process was guided by an independent editor, and the authors have also no other competing interests to declare.

**Disclaimer.** Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Acknowledgements.** This work utilized data from the OzFlux network, which is supported by the Australian Terrestrial Ecosystem Research Network (TERN; http://www.tern.org.au, last access: 8 March 2021). Lina Teckentrup, Martin G. De Kauwe, and Andrew J. Pitman acknowledge support from the Australian Research Council (ARC) Centre of Excellence for Climate Extremes (CE170100023). Martin G. De Kauwe and Andrew J. Pitman acknowledge support from the ARC Discovery Grant (DP190101823). Martin G. De Kauwe was also supported from the NSW Research Attraction and Acceleration Program. Sebastian Lienert acknowledges funding from the Swiss National Science Foundation (grant no. 172476) and from the European Union's Horizon 2020 research and innovation programme under grant agreement no. 821003 (project 4C, Climate-Carbon Interactions in the Current Century). ORNL is managed by UT-Battelle, LLC, for the DOE under contract DE-AC05-1008 00R22725. Emilie Joetzjer would like to thank the European Union's Horizon 2020 research and innovation programme with the CRESCENDO

project under the grant agreement no. 641816 and the H2020 CONSTRAIN under the grant agreement no. 820829. JSBACH simulations used resources of the Deutsches Klimarechenzentrum (DKRZ) under project ID bm0891. The CESM project is supported primarily by the National Science Foundation (NSF). This material is based upon work supported by the National Center for Atmospheric Research, which is a major facility sponsored by the NSF under cooperative agreement 1852977. Computing and data storage resources, including the Cheyenne supercomputer (https://doi.org/10.5065/D6RX99HX), were provided by the Computational and Information Systems Laboratory (CISL) at NCAR. Danica Lombardozzi was supported by the US Department of Agriculture NIFA award 2015-67003-23485. Sönke Zaehle was supported by the European Union's Horizon 2020 research and innovation programme under the grant agreement no. 821003 (4C project). Atul K. Jain acknowledges his funding support by the US Department of Energy (grant no. DE-SC0016323). Daniel S. Goll benefited from support from the Agence Nationale de la Recherche (ANR) grant ANR-16-CONV-0003 (CLAND). We thank Vladislav Bastrikov and Andrew J. Wiltshire for providing model output as part of the TRENDY v8 ensemble.

**Review statement.** This paper was edited by Ben Bond-Lamberty and reviewed by Ben Bond-Lamberty and one anonymous referee.

## References

Ahlström, A., Raupach, M. R., Schurgers, G., Smith, B., Arneth, A., Jung, M., Reichstein, M., Canadell, J. G., Friedlingstein, P., Jain, A. K., Kato, E., Poulter, B., Sitch, S., Stocker, B. D., Viovy, N., Wang, Y. P., Wiltshire, A., Zaehle, S., and Zeng, N.: The dominant role of semi-arid ecosystems in the trend and variability of the land CO₂ sink, Science, 348, 895–899, https://doi.org/10.1126/science.aaa1668, 2015.

Akagi, S. K., Yokelson, R. J., Wiedinmyer, C., Alvarado, M. J., Reid, J. S., Karl, T., Crounse, J. D., and Wennberg, P. O.: Emission factors for open and domestic biomass burning for use in atmospheric models, Atmos. Chem. Phys., 11, 4039–4072, https://doi.org/10.5194/acp-11-4039-2011, 2011.

Andreae, M. O. and Merlet, P.: Emission of trace gases and aerosols from biomass burning, Global Biogeochem. Cy., 15, 955–966, https://doi.org/10.1029/2000GB001382, https://doi.org/10.1029/2000GB001382, 2001.

Arora, V. K., Katavouta, A., Williams, R. G., Jones, C. D., Brovkin, V., Friedlingstein, P., Schwinger, J., Bopp, L., Boucher, O., Cadule, P., Chamberlain, M. A., Christian, J. R., Delire, C., Fisher, R. A., Hajima, T., Ilyina, T., Joetzjer, E., Kawamiya, M., Koven, C. D., Krasting, J. P., Law, R. M., Lawrence, D. M., Lenton, A., Lindsay, K., Pongratz, J., Raddatz, T., Séférrian, R., Tachiiri, K., Tjiputra, J. F., Wiltshire, A., Wu, T., and Ziehn, T.: Carbon-concentration and carbon-climate feedbacks in CMIP6 models and their comparison to CMIP5 models, Biogeosciences, 17, 4173–4222, https://doi.org/10.5194/bg-17-4173-2020, 2020.

Ballantyne, A. P., Alden, C. B., Miller, J. B., Tans, P. P., and White, J. W. C.: Increase in observed net carbon dioxide uptake by land and oceans during the past 50 years, Nature, 488, 70–72, https://doi.org/10.1038/nature11299, 2012.
Bastos, A., Running, S. W., Gouveia, C., and Trigo, R. M.: The global NPP dependence on ENSO: La Niña and the extraordinary year of 2011, J. Geophys. Res.-Biogeo., 118, 1247–1255, https://doi.org/10.1002/jgrg.20100, https://doi.org/10.1002/jgrg.20100, 2013.

Bastos, A., O'Sullivan, M., Ciais, P., Makowski, D., Sitch, S., Friedlingstein, P., Chevallier, F., Rödenbeck, C., Pongratz, J., Luijkx, I. T., Patra, P. K., Peylin, P., Canadell, J. G., Lauerwald, R., Li, W., Smith, N. E., Peters, W., Goll, D. S., Jain, A., Kato, E., Lienert, S., Lombardozzi, D. L., Haverd, V., Nabel, J. E. M. S., Poulter, B., Tian, H., Walker, A. P., and Zaehle, S.: Sources of Uncertainty in Regional and Global Terrestrial CO₂ Exchange Estimates, Global Biogeochem. Cy., 34, e2019GB006393, https://doi.org/10.1029/2019GB006393, 2020.

Beadle, N. C. W.: Soil Phosphate and Its Role in Molding Segments of the Australian Flora and Vegetation, with Special Reference to Xeromorphy and Sclerophylly, Ecology, 47, 992–1007, https://doi.org/10.2307/1935647, 1966.

Bloomfield, K. J., Prentice, I. C., Cernusak, L. A., Eamus, D., Medlyn, B. E., Rumman, R., Wright, I. J., Boer, M. M., Cale, P., Cleverly, J., Egerton, J. J. G., Ellsworth, D. S., Evans, B. J., Hayes, L. S., Hutchinson, M. F., Liddell, M. J., Macfarlane, C., Meyer, W. S., Togashi, H. F., Wardlaw, T., Zhu, L., and Atkin, O. K.: The validity of optimal leaf traits modelled on environmental conditions, New Phytol., 221, 1409–1423, https://doi.org/10.1111/nph.15495, 2019.

Box, E. O.: Plant functional types and climate at the global scale, J. Veg. Sci., 7, 309–320, https://doi.org/10.2307/3236274, 1996.

Brown, D., Taylor, J., and Bell, M.: The demography of desert Australia, Rangeland J., 30, 29–43, https://doi.org/10.1071/RJ07043, 2008.

Bureau of Meteorology, Commonwealth of Australia: Indian Ocean influences on Australian climate, available at: http://www.bom.gov.au/climate/iod/#tabs=Pacific-Ocean-interaction (last access: 25 August 2020), 2016.

Carvalhais, N., Forkel, M., Khomik, M., Bellarby, J., Jung, M., Migliavacca, M., Muu, M., Saatchi, S., Santoro, M., Thurner, M., Weber, U., Ahrens, B., Beer, C., Cescatti, A., Randerson, J., and Reichstein, M.: Global covariation of carbon turnover times with climate in terrestrial ecosystems, Nature, 514, 213–217, https://doi.org/10.1038/nature13731, 2014.

Cernusak, L. A., Hutley, L. B., Beringer, J., Holtum, J. A. M., and Turner, B. L.: Photosynthetic physiology of eucalypts along a sub-continental rainfall gradient in northern Australia, Agr. Forest Meteorol., 151, 1462–1470, https://doi.org/10.1016/j.agrformet.2011.01.006, savanna Patterns of Energy and Carbon Integrated Across the Landscape (SPECIAL), 2011.

Chapman, A.: Numbers of Living Species in Australia and the World, 1st Edn., Report for the Australian Biological Resources Study, Canberra, 2006.

Cleverly, J., Eamus, D., Luo, Q., Restrepo-Coupe, N., Kljun, N., Ma, X., Ewenz, C., Li, L., Yu, Q., and Huete, A.: The importance of interacting climate modes on Australia's contribution to global carbon cycle extremes, Sci. Rep., 6, 23113, https://doi.org/10.1038/srep23113, 2016.

Climate Prediction Center Internet Team: Historical El Niño / La Niña episodes (1950–present), available at:

https://origin.cpc.ncep.noaa.gov/products/analysis_monitoring/ensostuff/ONI_v5.php, last access: 25 August 2021.

De Kauwe, M. G., Medlyn, B. E., Ukkola, A. M., Mu, M., Sabot, M. E. B., Pitman, A. J., Meir, P., Cernusak, L. A., Rifai, S. W., Choat, B., Tissue, D. T., Blackman, C. J., Li, X., Roderick, M., and Briggs, P. R.: Identifying areas at risk of drought-induced tree mortality across South-Eastern Australia, Glob. Change Biol., 26, 5716–5733, https://doi.org/10.1111/gcb.15215, 2020.

Deser, C., Alexander, M. A., Xie, S.-P., and Phillips, A. S.: Sea Surface Temperature Variability: Patterns and Mechanisms, Ann. Rev. Mar. Sci., 2, 115–143, 21141660, https://doi.org/10.1146/annurev-marine-120408-151453, 2010.

DiMiceli, C. M., Carroll, M. L., Sohlberg, R. A., Huang, C., Hansen, M. C., and Townshend, J. R. G.: Annual global automated MODIS vegetation continuous fields (MOD44B) at 250 m spatial resolution for data years beginning day 65, 2000–2014, collection 5 percent tree cover, version 6, University of Maryland, College Park, MD, USA, available at: https://lpdaac.usgs.gov/products/mod44bv006/ (last access: 19 January 2021), 2017.

Donohue, R. J., McVicar, T. R., and Roderick, M. L.: Climate-related trends in Australian vegetation cover as inferred from satellite observations, 1981–2006, Glob. Change Biol., 15, 1025–1039, https://doi.org/10.1111/j.1365-2486.2008.01746.x, 2009.

Donohue, R. J., Roderick, M., McVicar, T. R., and Farquhar, G. D.: Impact of CO₂ fertilization on maximum foliage cover across the globe's warm, arid environments, Geophys. Res. Lett., 40, 3031–3035, https://doi.org/10.1002/grl.50563, 2013.

Drake, J. E., Aspinwall, M. J., Pfautsch, S., Rymer, P. D., Reich, P. B., Smith, R. A., Crous, K. Y., Tissue, D. T., Ghannoum, O., and Tjoelker, M.: The capacity to cope with climate warming declines from temperate to tropical latitudes in two widely distributed Eucalyptus species, Glob. Change Biol., 21, 459–472, https://doi.org/10.1111/gcb.12729, 2015.

Du, E., Terrer, C., Pellegrini, A. F. A., Ahlström, A., van Lissa, C. J., Zhao, X., Xia, N., Wu, X., and Jackson, R. B.: Global patterns of terrestrial nitrogen and phosphorus limitation, Nat. Geosci., 13, 221–226, https://doi.org/10.1038/s41561-019-0530-4, 2020.

ECMWF: CAMS Global Fire Assimilation System, available at: https://apps.ecmwf.int/datasets/data/cams-gfas/, last access: 6 March 2021.

Ellsworth, D., Anderson, I., Crous, K., Cooke, J., Drake, J., Gherlenda, A., Gimeno, T., Macdonald, C., Medlyn, B., Powell, J., Tjoelker, M., and Reich, P.: Elevated CO₂ does not increase eucalypt forest productivity on a low-phosphorus soil, Nat. Clim. Change, 7, 279–282, https://doi.org/10.1038/nclimate3235, 2017.

Fisher, R. A., Muszala, S., Verteinstein, M., Lawrence, P., Xu, C., McDowell, N. G., Knox, R. G., Koven, C., Holm, J., Rogers, B. M., Spessa, A., Lawrence, D., and Bonan, G.: Taking off the training wheels: the properties of a dynamic vegetation model without climate envelopes, CLM4.5(ED), Geosci. Model Dev., 8, 3593–3619, https://doi.org/10.5194/gmd-8-3593-2015, 2015.

Fleischer, K., Rammig, A., De Kauwe, M. G., Walker, A. P., Domingues, T. F., Fuchslueger, L., Garcia, S., Goll, D. S., Grandis, A., Jiang, M., Haverd, V., Hofhansl, F., Holm, J. A., Kruijt, B., Leung, F., Medlyn, B. E., Mercado, L. M., Norby, R. J., Pak, B., von Randow, C., Quesada, C. A., Schaap, K. J., Valverde-Barrantes, O. J., Wang, Y.-P., Yang, X., Zaehle, S., Zhu, Q., and Lapola, D. M.: Amazon forest response to CO₂ fertilization de-
pendent on plant phosphorus acquisition, Nat. Geosci., 12, 736-741, https://doi.org/10.1038/s41561-019-0404-9, 2019.

Friedlingstein, P., Jones, M. W., O'Sullivan, M., Andrew, R. M., Hauck, J., Peters, G. P., Peters, W., Pongratz, J., Sitch, S., Le Quéré, C., Bakker, D. C. E., Canadell, J. G., Ciais, P., Jackson, R. B., Anthoni, P., Barbero, L., Bastos, A., Bastrikov, V., Becker, M., Bopp, L., Buitenhuis, E., Chandra, N., Chevallier, F., Chini, L. P., Currie, K. I., Feely, R. A., Gehlen, M., Gilfillan, D., Gkritzalis, T., Goll, D. S., Gruber, N., Gutekunst, S., Harris, I., Haverd, V., Houghton, R. A., Hurtt, G., Ilyina, T., Jain, A. K., Joetzjer, E., Kaplan, J. O., Kato, E., Klein Goldewijk, K., Korsbakken, J. I., Landschützer, P., Lauvset, S. K., Lefèvre, N., Lenton, A., Lienert, S., Lombardozzi, D., Marland, G., McGuire, P. C., Melton, J. R., Metzl, N., Munro, D. R., Nabel, J. E. M. S., Nakaoka, S.-I., Neill, C., Omar, A. M., Ono, T., Peregon, A., Pierrot, D., Poulter, B., Rehder, G., Resplandy, L., Robertson, E., Rödenbeck, C., Séférían, R., Schwinger, J., Smith, N., Tans, P. P., Tian, H., Tilbrook, B., Tubiello, F. N., van der Werf, G. R., Wiltshire, A. J., and Zaehle, S.: Global Carbon Budget 2019, Earth Syst. Sci. Data, 11, 1783-1838, https://doi.org/10.5194/essd-11-1783-2019, 2019.

Friend, A. D., Lucht, W., Rademacher, T. T., Keribin, R., Betts, R., Cadule, P., Ciais, P., Clark, D. B., Dankers, R., Falloon, P. D., Ito, A., Kahana, R., Kleidon, A., Lomas, M. R., Nishina, K., Ostberg, S., Pavlick, R., Peylin, P., Schaphoff, S., Vuichard, N., Warszawski, L., Wiltshire, A., and Woodward, F. I.: Carbon residence time dominates uncertainty in terrestrial vegetation responses to future climate and atmospheric CO₂, P. Natl. Acad. Sci. USA, 111, 3280-3285, https://doi.org/10.1073/pnas.1222477110, 2014.

Hantson, S., Kelley, D. I., Arneth, A., Harrison, S. P., Archibald, S., Bachelet, D., Forrest, M., Hickler, T., Lasslop, G., Li, F., Mangeon, S., Melton, J. R., Nieradzik, L., Rabin, S. S., Prentice, I. C., Sheehan, T., Sitch, S., Teckentrup, L., Voulgarakis, A., and Yue, C.: Quantitative assessment of fire and vegetation properties in simulations with fire-enabled vegetation models from the Fire Model Intercomparison Project, Geosci. Model Dev., 13, 3299-3318, https://doi.org/10.5194/gmd-13-3299-2020, 2020.

Harris, I.: CRU TS v4.03: Climatic Research Unit (CRU) Time-Series (TS) version 4.03 of high-resolution gridded data of month-by-month variation in climate (Jan. 1901-Dec. 2018), Centre for Environmental Data Analysis (CEDA), https://doi.org/10.5285/10d3e3640f004c578403419aac167d82, 2019a.

Harris, I.: CRU JRA: Collection of CRU JRA forcing datasets of gridded land surface blend of Climatic Research Unit (CRU) and Japanese reanalysis (JRA) data, Centre for Environmental Data Analysis (CEDA), 2019b.

Harris, S. and Lucas, C.: Understanding the variability of Australian fire weather between 1973 and 2017, PLOS ONE, 14, e0222328, https://doi.org/10.1371/journal.pone.0222328, 2019.

Haverd, V., Raupach, M. R., Briggs, P. R., Canadell, J. G., Isaac, P., Pickett-Heaps, C., Roxburgh, S. H., van Gorsel, E., Viscarra Rossel, R. A., and Wang, Z.: Multiple observation types reduce uncertainty in Australia's terrestrial carbon and water cycles, Biogeosciences, 10, 2011-2040, https://doi.org/10.5194/bg-10-2011-2013, 2013a.

Haverd, V., Raupach, M. R., Briggs, P. R., J. G. Canadell., Davis, S. J., Law, R. M., Meyer, C. P., Peters, G. P., Pickett-Heaps, C., and

Sherman, B.: The Australian terrestrial carbon budget, Biogeosciences, 10, 851-869, https://doi.org/10.5194/bg-10-851-2013, 2013b.

Haverd, V., Smith, B., Raupach, M., Briggs, P., Nieradzik, L., Beringer, J., Hutley, L., Trudinger, C. M., and Cleverly, J.: Coupling carbon allocation with leaf and root phenology predicts tree-grass partitioning along a savanna rainfall gradient, Biogeosciences, 13, 761-779, https://doi.org/10.5194/bg-13-761-2016, 2016a.

Haverd, V., Smith, B., and Trudinger, C.: Process contributions of Australian ecosystems to interannual variations in the carbon cycle, Environ. Res. Lett., 11, 054013, https://doi.org/10.1088/1748-9326/11/5/054013, 2016b.

Haverd, V., Ahlström, A., Smith, B., and Canadell, J. G.: Carbon cycle responses of semi-arid ecosystems to positive asymmetry in rainfall, Glob. Change Biol., 23, 793-800, https://doi.org/10.1111/gcb.13412, 2017.

Haxeltine, A., Prentice, I. C., and Creswell, I. D.: A coupled carbon and water flux model to predict vegetation structure, J. Veg. Sci., 7, 651-666, https://doi.org/10.2307/3236377, 1996.

Hines, K. M., Bromwich, D. H., and Marshall, G. J.: Artificial Surface Pressure Trends in the NCEP-NCAR Reanalysis over the Southern Ocean and Antarctica, J. Climate, 13, 3940-3952, https://doi.org/10.1175/1520-0442(2000)013<3940:ASPTIT>2.0.CO;2, 2000.

Huntzinger, D. N., Michalak, A. M., Schwalm, C., Ciais, P., King, A. W., Fang, Y., Schaefer, K., Wei, Y., Cook, R. B., Fisher, J. B., Hayes, D., Huang, M., Ito, A., Jain, A. K., Lei, H., Lu, C., Maigan, F., Mao, J., Parazoo, N., Peng, S., Poulter, B., Ricciuto, D., Shi, X., Tian, H., Wang, W., Zeng, N., and Zhao, F.: Uncertainty in the response of terrestrial carbon sink to environmental drivers undermines carbon-climate feedback predictions, Sci. Rep., 7, 4765, https://doi.org/10.1038/s41598-017-03818-2, 2017.

Hurtt, G. C., Chini, L., Sahajpal, R., Frolking, S., Bodirsky, B. L., Calvin, K., Doelman, J. C., Fisk, J., Fujimori, S., Klein Goldewijk, K., Hasegawa, T., Havlik, P., Heinimann, A., Humpenöder, F., Jungclaus, J., Kaplan, J. O., Kennedy, J., Krisztin, T., Lawrence, D., Lawrence, P., Ma, L., Mertz, O., Pongratz, J., Popp, A., Poulter, B., Riahi, K., Shevliakova, E., Stehfest, E., Thornton, P., Tubiello, F. N., van Vuuren, D. P., and Zhang, X.: Harmonization of global land use change and management for the period 850-2100 (LUH2) for CMIP6, Geosci. Model Dev., 13, 5425-5464, https://doi.org/10.5194/gmd-13-5425-2020, 2020.

Hutchinson, M. F., McIntyre, S., Hobbs, R. J., Stein, J. L., Garnett, S., and Kinloch, J.: Integrating a global agro-climatic classification with bioregional boundaries in Australia, Global Ecol. Biogeogr., 14, 197-212, https://doi.org/10.1111/j.1466-822X.2005.00154.x, 2005.

IPCC: Managing the Risks of Extreme Events and Disasters to Advance Climate Change Adaptation: Special Report of the Intergovernmental Panel on Climate Change, Cambridge University Press, The Edinburgh Building, Shaftesbury Road, Cambridge CB2 8RU ENGLAND, 2012.

Jiang, M., Medlyn, B. E., Drake, J. E., Duursma, R. A., Anderson, I. C., Barton, C. V. M., Boer, M. M., Carrillo, Y., Castañeda-Gómez, L., Collins, L., Crous, K. Y., De Kauwe, M. G., dos Santos, B. M., Emmerson, K. M., Facey, S. L., Gherlenda, A. N., Gimeno, T. E., Hasegawa, S., Johnson, S. N., Kännaste, A., Mac-
donald, C. A., Mahmud, K., Moore, B. D., Nazaries, L., Neilson, E. H. J., Nielsen, U. N., Niinemets, Ü., Noh, N. J., Ochoa-Hueso, R., Pathare, V. S., Pendall, E., Pihblad, J., Piñeiro, J., Powell, J. R., Power, S. A., Reich, P. B., Renchon, A. A., Riegler, M., Rinnan, R., Rymer, P. D., Salomón, R. L., Singh, B. K., Smith, B., Tjoelker, M. G., Walker, J. K. M., Wujeska-Klause, A., Yang, J., Zaehle, S., and Ellsworth, D. S.: The fate of carbon in a mature forest under carbon dioxide enrichment, Nature, 580, 227–231, https://doi.org/10.1038/s41586-020-2128-9, 2020.

Joos, F. and Spahni, R.: Rates of change in natural and anthropogenic radiative forcing over the past 20,000 years, P. Natl. Acad. Sci. USA, 105, 1425–1430, https://doi.org/10.1073/pnas.0707386105, 2008.

Kaiser, J. W., Heil, A., Andreae, M. O., Benedetti, A., Chubarova, N., Jones, L., Morcrette, J.-J., Razinger, M., Schultz, M. G., Suttie, M., and van der Werf, G. R.: Biomass burning emissions estimated with a global fire assimilation system based on observed fire radiative power, Biogeosciences, 9, 527–554, https://doi.org/10.5194/bg-9-527-2012, 2012.

Keeling, R. and Manning, A.: Studies of Recent Changes in Atmospheric O₂ Content, Treatise on Geochemistry: Second Edition, 5, 385–404, https://doi.org/10.1016/B978-0-08-095975-7.00420-4, 2014.

Keenan, T. and Williams, C.: The Terrestrial Carbon Sink, Ann. Rev. Env. Resour., 43, 219–243, https://doi.org/10.1146/annurev-environ-102017-030204, 2018.

Kelley, D. and Harrison, S.: Enhanced Australian carbon sink despite increased wildfire during the 21st century, Environ. Res. Lett., 9, 104015, https://doi.org/10.1088/1748-9326/9/10/104015, 2014.

Kelley, D. I., Harrison, S. P., and Prentice, I. C.: Improved simulation of fire-vegetation interactions in the Land surface Processes and eXchanges dynamic global vegetation model (LPX-Mv1), Geosci. Model Dev., 7, 2411–2433, https://doi.org/10.5194/gmd-7-2411-2014, 2014.

Lambers, H., Finnegan, P. M., Jost, R., Plaxton, W. C., Shane, M. W., and Stitt, M.: Phosphorus nutrition in Proteaceae and beyond, Nat. Plants, 1, 15109, https://doi.org/10.1038/nplants.2015.109, 2015.

Le Quéré, C., Andrew, R. M., Friedlingstein, P., Sitch, S., Hauck, J., Pongratz, J., Pickers, P. A., Korsbakken, J. I., Peters, G. P., Canadell, J. G., Arneth, A., Arora, V. K., Barbero, L., Bastos, A., Bopp, L., Chevallier, F., Chini, L. P., Ciais, P., Doney, S. C., Gkritzalis, T., Goll, D. S., Harris, I., Haverd, V., Hoffman, F. M., Hoppema, M., Houghton, R. A., Hurtt, G., Ilyina, T., Jain, A. K., Johannessen, T., Jones, C. D., Kato, E., Keeling, R. F., Goldewijk, K. K., Landschützer, P., Lefèvre, N., Lienert, S., Liu, Z., Lombardozzi, D., Metzl, N., Munro, D. R., Nabel, J. E. M. S., Nakaoka, S., Neill, C., Olsen, A., Ono, T., Patra, P., Peregon, A., Peters, W., Peylin, P., Pfeil, B., Pierrot, D., Poulter, B., Rehder, G., Resplandy, L., Robertson, E., Rocher, M., Rödenbeck, C., Schuster, U., Schwinger, J., Séférion, R., Skjelvan, I., Steinhoff, T., Sutton, A., Tans, P. P., Tian, H., Tilbrook, B., Tubiello, F. N., van der Laan-Luijkx, I. T., van der Werf, G. R., Viovy, N., Walker, A. P., Wiltshire, A. J., Wright, R., Zaehle, S., and Zheng, B.: Global Carbon Budget 2018, Earth Syst. Sci. Data, 10, 2141–2194, https://doi.org/10.5194/essd-10-2141-2018, 2018.

Li, F., Val Martin, M., Andreae, M. O., Arneth, A., Hantson, S., Kaiser, J. W., Lasslop, G., Yue, C., Bachelet, D., Forrest, M.,

Kluzek, E., Liu, X., Mangeon, S., Melton, J. R., Ward, D. S., Darmenov, A., Hickler, T., Ichoku, C., Magi, B. I., Sitch, S., van der Werf, G. R., Wiedinmyer, C., and Rabin, S. S.: Historical (1700–2012) global multi-model estimates of the fire emissions from the Fire Modeling Intercomparison Project (FireMIP), Atmos. Chem. Phys., 19, 12545–12567, https://doi.org/10.5194/acp-19-12545-2019, 2019.

Li, X. and Xiao, J.: A Global, 0.05-Degree Product of Solar-Induced Chlorophyll Fluorescence Derived from OCO-2, MODIS, and Reanalysis Data, Remote Sens., 11, 517, https://doi.org/10.3390/rs11050517, 2019 (data available at: https://globalecology.unh.edu/data/GOSIF-GPP.html, last access: 20 April 2020).

Li, X., Blackman, C. J., Choat, B., Duursma, R. A., Rymer, P. D., Medlyn, B. E., and Tissue, D. T.: Tree hydraulic traits are coordinated and strongly linked to climate-of-origin across a rainfall gradient, Plant Cell Environ., 41, 646–660, https://doi.org/10.1111/pce.13129, 2018.

Liu, Y. Y., van Dijk, A. I. J. M., de Jeu, R. A. M., Canadell, J. G., McCabe, M. F., Evans, J. P., and Wang, G.: Recent reversal in loss of global terrestrial biomass, Nat. Clim. Change, 5, 470–474, https://doi.org/10.1038/nclimate2581, 2015 (data available at: http://wald.anu.edu.au/global-biomass/, last access: 13 August 2020).

Luo, Y.: Terrestrial Carbon-Cycle Feedback to Climate Warming, Annu. Rev. Ecol. Evol. S., 38, 683–712, https://doi.org/10.1146/annurev.ecolsys.38.091206.095808, 2007.

Luo, Y., White, L. W., Canadell, J. G., DeLucia, E. H., Ellsworth, D. S., Finzi, A., Lichter, J., and Schlesinger, W. H.: Sustainability of terrestrial carbon sequestration: A case study in Duke Forest with inversion approach, Global Biogeochem. Cy., 17, 1021, https://doi.org/10.1029/2002GB001923, 2003.

Ma, L., Hurtt, G. C., Chini, L. P., Sahajpal, R., Pongratz, J., Frolking, S., Stehfest, E., Klein Goldewijk, K., O'Leary, D., and Doelman, J. C.: Global rules for translating land-use change (LUH2) to land-cover change for CMIP6 using GLM2, Geosci. Model Dev., 13, 3203–3220, https://doi.org/10.5194/gmd-13-3203-2020, 2020.

Ma, X., Huete, A., Cleverly, J., Eamus, D., Chevallier, F., Joanna, J., Poulter, B., Zhang, Y., Guanter, L., Meyer, W., Xie, Z., and Ponce-Campos, G.: Drought rapidly diminishes the large net CO₂ uptake in 2011 over semi-arid Australia, Sci. Rep., 6, 37747, https://doi.org/10.1038/srep37747, 2016.

Marquer, L., Dallmeyer, A., Poska, A., Pongratz, J., Smith, B., and Gaillard, M.-J.: Modeling past human-induced vegetation change is a challenge – the case of Europe, Past Global Change Magazine, 26, 12–13, https://doi.org/10.22498/pages.26.1.12, 2018.

Medlyn, B. E., De Kauwe, M. G., Zaehle, S., Walker, A. P., Duursma, R. A., Luus, K., Mishurov, M., Pak, B., Smith, B., Wang, Y.-P., Yang, X., Crous, K. Y., Drake, J. E., Gimeno, T. E., Donald, C. A., Norby, R. J., Power, S. A., Tjoelker, M. G., and Ellsworth, D. S.: Using models to guide field experiments: a priori predictions for the CO₂ response of a nutrient- and water-limited native Eucalypt woodland, Glob. Change Biol., 22, 2834–2851, https://doi.org/10.1111/gcb.13268, 2016.

Meiyappan, P., Jain, A. K., and House, J. I.: Increased influence of nitrogen limitation on CO₂ emissions from future land use
and land use change, Global Biogeochem. Cy., 29, 1524–1548,
https://doi.org/10.1002/2015GB005086, 2015.

Meyerholt, J., Sickel, K., and Zaehle, S.: Ensemble projections elucidate effects of uncertainty in terrestrial nitrogen limitation on future carbon uptake, Glob. Change Biol., 26, 3978–3996, https://doi.org/10.1111/gcb.15114, 2020.

Mitchell, P. J., O'Grady, A. P., Hayes, K. R., and Pinkard, E. A.: Exposure of trees to drought-induced die-off is defined by a common climatic threshold across different vegetation types, Ecol. Evol., 4, 1088–1101, https://doi.org/10.1002/ece3.1008, 2014.

OzFlux: Australian and New Zealand Flux Research and Monitoring, available at: http://www.ozflux.org.au/, last access: 26 April 2017.

Pan, X., Ichoku, C., Chin, M., Bian, H., Darmenov, A., Colarco, P., Ellison, L., Kucsera, T., da Silva, A., Wang, J., Oda, T., and Cui, G.: Six global biomass burning emission datasets: intercomparison and application in one global aerosol model, Atmos. Chem. Phys., 20, 969–994, https://doi.org/10.5194/acp-20-969-2020, 2020.

Piao, S., Sitch, S., Ciais, P., Friedlingstein, P., Peylin, P., Wang, X., Ahlström, A., Anav, A., Canadell, J. G., Cong, N., Huntingford, C., Jung, M., Levis, S., Levy, P. E., Li, J., Lin, X., Lomas, M. R., Lu, M., Luo, Y., Ma, Y., Myneni, R. B., Poulter, B., Sun, Z., Wang, T., Viovy, N., Zaehle, S., and Zeng, N.: Evaluation of terrestrial carbon cycle models for their response to climate variability and to CO₂ trends, Glob. Change Biol., 19, 2117–2132, https://doi.org/10.1111/gcb.12187, 2013.

Pongratz, J., D., H., Don, A., Erb, K.-H., Fuchs, R., Herold, M., Jones, C., Kuemmerle, T., Luyssaert, S., Meyfroidt, P., and Naudts, K.: Models meet data: Challenges and opportunities in implementing land management in Earth system models, Glob. Change Biol., 24, 1470–1487, https://doi.org/10.1111/gcb.13988, 2018.

Poulter, B., Frank, D., Ciais, P., Myneni, R., Andela, N., Bi, J., Broquet, G., Canadell, J., Chevallier, F., Liu, Y., Running, S., Sitch, S., and van der Werf, G.: Contribution of semi-arid ecosystems to interannual variability of the global carbon cycle, Nature, 509, 600–603, https://doi.org/10.1038/nature13376, 2014.

Pugh, T. A. M., Lindeskog, M., Smith, B., Poulter, B., Arneth, A., Haverd, V., and Calle, L.: Role of forest regrowth in global carbon sink dynamics, P. Natl. Acad. Sci. USA, 116, 4382–4387, https://doi.org/10.1073/pnas.1810512116, 2019.

Pugh, T. A. M., Rademacher, T., Shafer, S. L., Steinkamp, J., Barichivich, J., Beckage, B., Haverd, V., Harper, A., Heinke, J., Nishina, K., Rammig, A., Sato, H., Arneth, A., Hantson, S., Hickler, T., Kautz, M., Quesada, B., Smith, B., and Thonicke, K.: Understanding the uncertainty in global forest carbon turnover, Biogeosciences, 17, 3961–3989, https://doi.org/10.5194/bg-17-3961-2020, 2020.

Raupach, M. R., Canadell, J. G., and Le Quéré, C.: Anthropogenic and biophysical contributions to increasing atmospheric CO₂ growth rate and airborne fraction, Biogeosciences, 5, 1601–1613, https://doi.org/10.5194/bg-5-1601-2008, 2008.

Rogers, A.: The use and misuse of $V_{c,max}$ in Earth System Models, Photosynth. Res., 119, 15–29, https://doi.org/10.1007/s11120-013-9818-1, 2013.

Roxburgh, S., Barrett, D., Berry, S., Carter, J., Davies, I., Gifford, R., Kirschbaum, M., McBeth, B., Noble, I., Parton, W., Raupach, M., and Roderick, M.: A critical overview of model estimates of

net primary productivity for the Australian continent, Funct Plant
Biol., 31, 1043–1059, https://doi.org/10.1071/FP04100, 2004.

Sakschewski, B., von Bloh, W., Drüke, M., Sörensson, A. A., Rus-
cica, R., Langerwisch, F., Billing, M., Bereswill, S., Hirota,
M., Oliveira, R. S., Heinke, J., and Thonicke, K.: Variable tree
rooting strategies are key for modelling the distribution, pro-
ductivity and evapotranspiration of tropical evergreen forests,
Biogeosciences, 18, 4091–4116, https://doi.org/10.5194/bg-18-
4091-2021, 2021.

Scott, D. W.: Multivariate density estimation: theory, practice, and visualization, Wiley, New York, 1992.

Smith, B., Wårlind, D., Arneth, A., Hickler, T., Leadley, P., Siltberg, J., and Zaehle, S.: Implications of incorporating N cycling and N limitations on primary production in an individual-based dynamic vegetation model, Biogeosciences, 11, 2027–2054, https://doi.org/10.5194/bg-11-2027-2014, 2014.

Tarin, T., Nolan, R. H., Eamus, D., and Cleverly, J.: Carbon and water fluxes in two adjacent Australian semi-arid ecosystems, Agr. Forest Meteorol., 281, 107853, https://doi.org/10.1016/j.agrformet.2019.107853, 2020.

Teckentrup, L.: Iteckentrup/TRENDY_v8_Australia: v1.0.0 (v1.0.0), Zenodo [code], https://doi.org/10.5281/zenodo.5570974, 2021.

Teckentrup, L., Harrison, S. P., Hantson, S., Heil, A., Melton, J. R., Forrest, M., Li, F., Yue, C., Arneth, A., Hickler, T., Sitch, S., and Lasslop, G.: Response of simulated burned area to historical changes in environmental and anthropogenic factors: a comparison of seven fire models, Biogeosciences, 16, 3883–3910, https://doi.org/10.5194/bg-16-3883-2019, 2019.

Thomas, R. Q., Zaehle, S., Templer, P. H., and Goodale, C. L.: Global patterns of nitrogen limitation: confronting two global biogeochemical models with observations, Glob. Change Biol., 19, 2986–2998, https://doi.org/10.1111/gcb.12281, 2013.

Tian, H., Yang, J., Lu, C., Xu, R., Canadell, J. G., Jackson, R. B., Ar-
neth, A., Chang, J., Chen, G., Ciais, P., Gerber, S., Ito, A., Huang,
Y., Joos, F., Lienert, S., Messina, P., Olin, S., Pan, S., Peng,
C., Saikawa, E., Thompson, R. L., Vuichard, N., Winiwarter,
W., Zaehle, S., Zhang, B., Zhang, K., and Zhu, Q.: The Global
N₂O Model Intercomparison Project, B. Am. Meteorol. Soc., 99,
1231–1251, https://doi.org/10.1175/BAMS-D-17-0212.1, 2018.

Togashi, H. F., Prentice, I. C., Evans, B. J., Forrester, D. I., Drake, P., Feikema, P., Brooksbank, K., Eamus, D., and Taylor, D.: Morphological and moisture availability controls of the leaf area-to-sapwood area ratio: analysis of measurements on Australian trees, Ecol. Evol., 5, 1263–1270, https://doi.org/10.1002/ece3.1344, 2015.

Trancoso, R., Larsen, J. R., McVicar, T. R., Phinn, S. R., and McAlpine, C. A.: CO₂-vegetation feedbacks and other climate changes implicated in reducing base flow, Geophys. Res. Lett., 44, 2310–2318, https://doi.org/10.1002/2017GL072759, 2017.

Ukkola, A. M., De Kauwe, M. G., Pitman, A. J., Best, M. J., Abramowitz, G., Haverd, V., Decker, M., and Haughton, N.: Land surface models systematically overestimate the intensity, duration and magnitude of seasonal-scale evaporative droughts, Environ. Res. Lett., 11, 104012, https://doi.org/10.1088/1748-9326/11/10/104012, 2016a.

Ukkola, A. M., Prentice, I. C., Keenan, T. F., van Dijk, A. I. J. M., Viney, N. R., Myneni, R. B., and Bi, J.: Reduced streamflow in water-stressed climates consistent with
CO₂ effects on vegetation, Nat. Clim. Change, 6, 75–78, https://doi.org/10.1038/nclimate2831, 2016b.

Ummenhofer, C. C., Sen Gupta, A., Briggs, P. R., England, M. H., McIntosh, P. C., Meyers, G. A., Pook, M. J., Raupach, M. R., and Risbey, J. S.: Indian and Pacific Ocean Influences on Southeast Australian Drought and Soil Moisture, J. Climate, 24, 1313–1336, https://doi.org/10.1175/2010JCLI3475.1, 2011.

van der Horst, S. V. J., Pitman, A. J., De Kauwe, M. G., Ukkola, A., Abramowitz, G., and Isaac, P.: How representative are FLUXNET measurements of surface fluxes during temperature extremes?, Biogeosciences, 16, 1829–1844, https://doi.org/10.5194/bg-16-1829-2019, 2019.

van der Werf, G. R., Randerson, J. T., Giglio, L., van Leeuwen, T. T., Chen, Y., Rogers, B. M., Mu, M., van Marle, M. J. E., Morton, D. C., Collatz, G. J., Yokelson, R. J., and Kasibhatla, P. S.: Global fire emissions estimates during 1997–2016, Earth Syst. Sci. Data, 9, 697–720, https://doi.org/10.5194/essd-9-697-2017, 2017 (data available at: https://globalfiredata.org/pages/data/, last access: 11 July 2020).

van Gorsel, E., Wolf, S., Cleverly, J., Isaac, P., Haverd, V., Ewenz, C., Arndt, S., Beringer, J., Resco de Dios, V., Evans, B. J., Griebel, A., Hutley, L. B., Keenan, T., Kljun, N., Macfarlane, C., Meyer, W. S., McHugh, I., Pendall, E., Prober, S. M., and Silberstein, R.: Carbon uptake and water use in woodlands and forests in southern Australia during an extreme heat wave event in the “Angry Summer” of 2012/2013, Biogeosciences, 13, 5947–5964, https://doi.org/10.5194/bg-13-5947-2016, 2016.

Walker, A. P., De Kauwe, M. G., Bastos, A., Belmecheri, S., Georgiou, K., Keeling, R., McMahon, S. M., Medlyn, B. E., Moore, D. J., Norby, R. J., Zaehle, S., Anderson-Teixeira, K. J., Battipaglia, G., Brienen, R. J., Cabugao, K. G., Cailleret, M., Campbell, E., Canadell, J., Ciais, P., Craig, M. E., Ellsworth, D., Farquhar, G., Fatichi, S., Fisher, J. B., Frank, D., Graven, H., Gu, L., Haverd, V., Heilman, K., Heimann, M., Hungate, B. A., Iversen, C. M., Joos, F., Jiang, M., Keenan, T. F., Knauer, J., Körner, C., Leshyk, V. O., Leuzinger, S., Liu, Y., MacBean, N., Malhi, Y., McVicar, T. R., Penuelas, J., Pongratz, J., Powell, A. S., Riutta, T., Sabot, M. E. B., Schleucher, J., Sitch, S., Smith, W. K., Sulman, B., Taylor, B., Terrer, C., Torn, M. S., Treseder, K. K., Trugman, A. T., Trumbore, S. E., van Mantgem, P. J., Voelker, S. L., Whelan, M. E., and Zuidema, P. A.: Integrating the evidence for a terrestrial carbon sink caused by increasing atmospheric CO₂, New Phytol., https://doi.org/10.1111/nph.16866, 2020.

Whitley, R., Beringer, J., Hutley, L. B., Abramowitz, G., De Kauwe, M. G., Duursma, R., Evans, B., Haverd, V., Li, L., Ryu, Y., Smith, B., Wang, Y.-P., Williams, M., and Yu, Q.: A model inter-comparison study to examine limiting factors in modelling Australian tropical savannas, Biogeosciences, 13, 3245–3265, https://doi.org/10.5194/bg-13-3245-2016, 2016.

Whitley, R., Beringer, J., Hutley, L. B., Abramowitz, G., De Kauwe, M. G., Evans, B., Haverd, V., Li, L., Moore, C., Ryu, Y., Scheiter, S., Schymanski, S. J., Smith, B., Wang, Y.-P., Williams, M., and Yu, Q.: Challenges and opportunities in land surface modelling of savanna ecosystems, Biogeosciences, 14, 4711–4732, https://doi.org/10.5194/bg-14-4711-2017, 2017.

Wild, A.: The phosphate content of Australian soils, Aust. J. Agr. Res., 9, 193–204, https://doi.org/10.1071/ar9580193, 1958.

Yang, J., Medlyn, B. E., De Kauwe, M. G., Duursma, R. A., Jiang, M., Kumarathunge, D., Crous, K. Y., Gimeno, T. E., Wujeska-Klause, A., and Ellsworth, D. S.: Low sensitivity of gross primary production to elevated CO₂ in a mature eucalypt woodland, Biogeosciences, 17, 265–279, https://doi.org/10.5194/bg-17-265-2020, 2020.

Zaehle, S.: Terrestrial nitrogen-carbon cycle interactions at the global scale, Philos. T. R. Soc. B, 368, 20130125, https://doi.org/10.1098/rstb.2013.0125, 2013.

Zhang, X., Wang, Y.-P., Peng, S., Rayner, P. J., Ciais, P., Silver, J. D., Piao, S., Zhu, Z., Lu, X., and Zheng, X.: Dominant regions and drivers of the variability of the global land carbon sink across timescales, Glob. Change Biol., 24, 3954–3968, https://doi.org/10.1111/gcb.14275, 2018.

Zhang, Y., Dannenberg, M. P., Hwang, T., and Song, C.: El Niño-Southern Oscillation-Induced Variability of Terrestrial Gross Primary Production During the Satellite Era, J. Geophys. Res.-Biogeo., 124, 2419–2431, https://doi.org/10.1029/2019JG005117, 2019.