Process contributions of Australian ecosystems to interannual variations in the carbon cycle

New evidence is emerging that semi-arid ecosystems dominate interannual variability (IAV) of the global carbon cycle, largely via fluctuating water availability associated with El Niño/Southern Oscillation. Recent evidence from global terrestrial biosphere modelling and satellite-based inversion of atmospheric CO2 point to a large role of Australian ecosystems in global carbon cycle variability, including a large contribution from Australia to the record land sink of 2011. However the specific mechanisms governing this variability, and their bioclimatic distribution within Australia, have not been identified. Here we provide a regional assessment, based on best available observational data, of IAV in the Australian terrestrial carbon cycle and the role of Australia in the record land sink anomaly of 2011. We find that IAV in Australian net carbon uptake is dominated by semi-arid ecosystems in the east of the continent, whereas the 2011 anomaly was more uniformly spread across most of the continent. Further, and in contrast to global modelling results suggesting that IAV in Australian net carbon uptake is amplified by lags between production and decomposition, we find that, at continental scale, annual variations in production are dampened by annual variations in decomposition, with both fluxes responding positively to precipitation anomalies.


Introduction
There is compelling new evidence of the important role of semi-arid ecosystems in the variability of the global carbon cycle. In a regional attribution of terrestrial biosphere model predictions of global net biome productivity (NBP), Ahlström et al (2015) demonstrated that, while tropical forests dominate the mean of global NBP , semi-arid ecosystems dominate both its interannual variability (IAV) and trend. The same study found that in semi-arid ecosystems, IAV in modelled NBP is dominated by the response of vegetation productivity to variations in water availability, strongly associated with the El Niño/Southern Oscillation (ENSO). ENSO explains more than 40% of global net primary production (NPP) variability, mainly driven by the response of Southern Hemisphere ecosystems, particularly in tropical and subtropical regions (Bastos et al 2013). The large role of semi-arid ecosystems in global carbon uptake anomalies is further exemplified by the significant contribution of Australia, the driest inhabited continent, to the record residual global land sink of 2011 of 1.5±0.9 Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence.
Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. Detmers et al (2015) also inferred a large Australian land sink in 2011 (0.79 PgC).
These earlier studies point out that 2011 was an exceptional year in terms of the regional climate and its effects on ecosystem functioning, but the specific climatic drivers and responding mechanisms have not been clearly identified. In this work we upscale carbon and water fluxes from 14 sites of the OzFlux network (Isaac 2014) using a biogeochemical land surface model further constrained by observations of biomass, soil carbon, streamflow and remotely sensed vegetation cover (Zhu et al 2013) (updated from Haverd et al 2013b). We use the results to quantify IAV in the Australian terrestrial carbon cycle, including the magnitude of the 2011 anomaly, and go on to attribute the IAV and its major component fluxes spatially and to responding processes.

Methods
Australian net ecosystem productivity (NEP)  and its components, R H and NPP, were derived using model-data synthesis (Trudinger et al 2016), updated from (Haverd et al 2013a), the only published assessments of Australian biospheric carbon balance which integrate multiple regional observation sources. Of these sources, remotely sensed vegetation cover and flux-tower observations of carbon and water exchange are particularly relevant constraints on interannual variations. Evaluation metrics (Trudinger et al 2016) reveal good prediction of OzFlux annual gross primary production (n=70, r 2 =0.63, normalised root mean squared error ( ) s y x , x 2 / NRMSE=0.4), with poor correlation but acceptable variance in NEP (n=70, r 2 =0.01, NRMSE=1.1).
The model-data synthesis is described briefly below: full details are provided by Haverd et al (2013a) and Trudinger et al (2016). The synthesis employs BIOS2, a fine-spatial-resolution (0.05°) offline modelling environment, including a modification of the CABLE biogeochemical land surface model (Wang et al 2010(Wang et al , 2011b incorporating the SLI soil model (Haverd and Cuntz 2010). BIOS2 parameters are constrained and predictions are evaluated using multiple observation sets from across the Australian continent, including streamflow from 416 gauged catchments, eddy flux data (CO 2 and H 2 O) from 14 OzFlux sites (Isaac 2014), litterfall data, and soil, litter and biomass carbon pools.
CABLE consists of five components (Wang et al 2011a): (1) the radiation module describes direct and diffuse radiation transfer and absorption by sunlit and shaded leaves; (2) the canopy micrometeorology module describes the surface roughness length, zeroplane displacement height, and aerodynamic conductance from the reference height to the air within canopy or to the soil surface; (3) the canopy module includes the coupled energy balance, transpiration, stomatal conductance and photosynthesis of sunlit and shaded leaves; (4) the soil module describes heat and water fluxes within soil and snow at their respective surfaces; and (5) the ecosystem carbon module accounts for the respiration of stem, root and soil organic carbon decomposition. In BIOS2, the default CABLE v1.4 soil and carbon modules were replaced respectively by the SLI soil model ( Parameter estimation is achieved using the Levenberg-Marquardt method, as implemented in the PEST software (Doherty et al 2010) to minimise residuals between predictions and observations. An ensemble of parameter sets that are consistent with the observations are generated using the null space Monte Carlo method (Doherty et al 2010). Variance of predictions based on this ensemble of parameter sets is equated with variance attributable to parameter equifinality in model predictions.
Vegetation cover is prescribed, using a monthly time series of fractional absorbed photosynthetically active radiation (fPAR) for January 1982 to December 2013 that was derived from the third generation . A monthly maximum composite was created from the original 15 day series, and the data were resampled from the original 0.0833°resolution (8 km) to 0.05°(5 km). NDVI values from 0.1 (bare ground) to 0.75 (full cover) were linearly rescaled between 0 and 1 to represent vegetation fractional cover. Total fPAR is partitioned into persistent (mainly woody) and recurrent (mainly grassy) vegetation components, following the methodology of Donohue et al (2009) and Lu et al (2003). This methodology takes advantage of low levels of seasonal change in LAI in woody vegetation, allowing seasonal variation in fPAR to be attributed principally to grassy vegetation. The remaining and relatively constant fPAR signal is attributed to woody vegetation. LAI for woody and grassy components are estimated by Beer's Law (e.g. Houldcroft et al 2009): where the vegetation type is either W (persistent or mainly woody) or G (recurrent or mainly grassy) and k is an extinction coefficient, set here to 0.5. In contrast to earlier BIOS2 simulations (Haverd et al 2013a), equation (2) accounts for the effect of shading of grass by woody vegetation. We partition flux variability amongst component fluxes (components may be spatial or process contributions) using the formalism of Ahlström et al (2015) (equation (3)) Here x jt is the flux anomaly (departure from a longterm trend) for the jth constituent flux at time t (in years), and X t is the total flux anomaly, with = å X x .
t t jt By this definition, the jth contribution to the total IAV, f j , is the average relative anomaly x jt /X t for the jth constituent, weighted with the absolute global anomaly |Xt|. The definition ensures that å = f 1, j j but allows individual contributions to fall outside the range (0, 1) if the total anomaly X t arises from partially cancelling contributions x jt from different constituents.
We focus on IAV in NEP, which is the difference between NPP and heterotrophic respiration (R H ), in the absence of disturbance. We do not account here for the component of IAV in NBP attributable to disturbance (land use and natural disturbance agents, mainly fire): the IAV of gross Australian fire emissions (based on the Global Fire Emissions Database (GFED3); van der Werf et al 2010) is less than 20% of IAV in Australian NEP (Haverd et al 2013b), while global modelling studies suggest that the net contribution of wildfires to semi-arid ecosystem IAV is small (Ahlström et al 2015). We confirm this by presenting updated gross annual fire emission anomalies (GFED4s) (udated from van der Werf et al 2010) alongside our estimated NEP anomalies.

Spatial attribution of IAV of continental net carbon uptake and its component fluxes
For the purpose of regional attribution, we use the bioclimatic classification of figure 2(i) which is an aggregation of classes from the agro-climatic classification of Hutchinson et al (2005) (table 2, figure 3). Australian biogeography strongly reflects rainfall pattern ( figure 2(ii)). The contributions (equation (1)) of the ecosystems of each region to IAV in continental NPP, R H and NEP are shown in figure 2(iii). The Sparsely Vegetated and Savanna regions dominate contributions to continental inter-annual variability in all three fluxes, together representing a 90% contribution to IAV in continental NEP. Notably, the Sparsely Vegetated inland region, which contributes 36% to IAV in continental NPP, has a much larger contribution of 46% to variability in continental NEP, a result we will explore further in section 6 when we quantify process contributions to IAV in NEP.
Figures 2(iv)-(ix) shows the spatial distribution of IAV in NPP, R H and NEP fluxes, and the spatial contributions to IAV in the continental fluxes. At the annual time-scale, NPP is most variable across the eastern part of the savanna belt (figure 2(iv)), a region of highly variable rainfall where vegetation is both productive and strongly responsive to water availability (Raupach et al 2013). Variability of heterotrophic respiration (figure 2(v)) is high along the eastern seaboard where high labile carbon stores coincide with high variability in soil moisture availability. The variability in NEP (figure 2(vi)) is lower overall than the constituent fluxes and with a spatial distribution reflecting elements of both components. The spatial contribution maps (figures 2(vii)-(ix)) contrast with the maps of IAV because of spatial differences in the timing of the variability. For example, significant IAV in C fluxes in the South-West of Western Australia contribute negligibly to continental IAV, suggesting that the variability in this region is not correlated with that of the continent as a whole, a pattern also evident for the precipitation driver (not shown), and related to the dominant influence of the Indian Ocean Dipole on rainfall variability of the South-West, while the ENSO is the key driver of rainfall and weather patterns over much of the continent (Risbey et al 2009).

Extremeness of Australian carbon cycle in 2011
The extremeness of the 2011 carbon-cycle anomaly is quantified in figure 3 by the z-score (number of standard deviations from the local 1982-2013 mean flux). The regionally averaged NPP anomaly was anomalously high (z-score>1) across all bioclimatic regions, and particularly high (>3) in the semi-arid regions. Heterotrophic respiration was similarly elevated across all regions, although less extreme than NPP in the semi-arid regions, leading to particularly elevated NEP (z-score>3) in these regions and for Australia as a whole. The spatial distributions (figures 3(ii)-(iv)) show the widespread extent of the anomaly.
The spatial pattern of flux anomalies depends not only on the local z-score, but on the local average flux,

Process contributions to IAV in net carbon uptake
Thus far our analysis demonstrates a dampening of interannual variations in continental NEP compared to NPP (figures 1, 2(iv), (vi)) because of temporal correlations between production (NPP) and decomposition (R H ). We now quantify the contributions (equation (3)) of the constituent fluxes (NPP and R H ) to IAV. Figure 4(i) shows the contributions of regional NPP and R H to continental IAV in NEP. A positive contribution indicates that the constituent flux is correlated with continental NEP, and a negative contribution that it is anti-correlated with continental NEP at the annual time scale. Arrows indicate the sign and magnitude of the relative regional contributions to continental NEP variability. Across all regions, and for the continent as a whole, we see that variations in NPP predominate, but are significantly compensated for by variations in R H in terms of their contribution to continental NEP variability. R H variability offsets the contribution of NPP variability by 40% continentally: regionally, the offset is largest in the Savanna region.
The spatial patterns of the NEP and R H contributions to IAV in NEP (figures 4(ii) and (iii)) reinforce the picture of a dominant contribution by NPP, significantly offset by R H across most of the continent. This is because both constituent fluxes respond positively to precipitation anomalies. Parts of the cool temperate regions (Tasmania and Southern Victoria) are exceptional, being characterised by IAV in the R H component that re-inforces the NPP component (i.e. positive NPP anomalies occur in the same years as negative R H anomalies).

Conclusion
Our results provide a regional assessment, based on best available observational data, of IAV in the Australian terrestrial carbon cycle and the role of Australia in the record land sink anomaly of 2011. Our results fall in line with Ahlström et al (2015), Poulter et al (2014), Detmers et al (2015) and Bastos et al (2013), confirming Australia's large role in the record land sink anomaly of 2011. Our results demonstrate that such variability is consistent with a direct physiological response of vegetation productivity to fluctuating water availability, with these interannual variations in productivity being partially offset by a largely correlated opposing effect of fluctuating water availability on decomposition.