Skillful seasonal prediction of key carbon cycle components: NPP and fire risk

We investigate the skill of the GloSea5 seasonal forecasting system for two carbon cycle processes, which are strong contributors to global CO2 variability: the impact of meteorological conditions on CO2 uptake by vegetation (characterised by net primary productivity, NPP), and on fire occurrences (characterised by fire risk indices). Current seasonal forecasts of global CO2 concentrations rely on the relationship with the El Niño–Southern Oscillation (ENSO), combined with estimated anthropogenic emissions. NPP and fire are key processes underlying that global CO2–ENSO relationship: In the tropics, during El Niño events, CO2 uptake by vegetation is reduced and fires occur more frequently, leading to higher global CO2 levels. Our study assesses the skill of these processes in the forecast model for the first time. We use the McArthur forest fire index, calculated from daily data from several meteorological variables. We also assess a simpler fire index, based solely on seasonal mean temperature and relative humidity, to test the need for additional complexity. For NPP, the skill is high in regions that respond strongly to ENSO, such as equatorial South America in boreal winter, and northeast Brazil in boreal summer. There is also skill in some regions without a strong ENSO response. The fire risk indices show significant skill across much of the tropics, including Indonesia, southern and eastern Africa, and parts of the Amazon. We relate this skill to the underlying meteorological variables, finding that fire risk in particular follows similar patterns to relative humidity. On the seasonal-mean timescale, the McArthur index offers no benefits over the simpler fire index: they show the same relationship to burnt area and response to ENSO, and the same levels of skill, in almost all cases. Our results highlight potentially useful prediction skill, as well as important limitations, for seasonal forecasts of land-surface impacts of climate variability.


Introduction
Atmospheric CO 2 levels are often considered to be a 'driver' of changes in the climate system, particularly when considering projected changes over many decades. However CO 2 concentration is itself a function of both anthropogenic emissions and exchanges with the land and ocean. It is these natural fluxes, and in particular the three factors of heterotrophic respiration, vegetation productivity and fire variation, that dominate the interannual variability of CO 2 (Zeng et al 2005). It has been shown that CO 2 concentrations can be skilfully predicted from features of the climate system, such as the El Niño-Southern Oscillation (ENSO) together with an estimate of anthropogenic emissions (Jones et al 2001, Jones and Cox 2005. However, seasonal forecasts based on global climate models could also be used to forecast components of CO 2 variability more directly. We aim here to investigate that possibility, for vegetation productivity and fire risk, in the context of the skillful ENSO-based CO 2 forecasts.
ENSO is the largest global mode of interannual variability and affects both the global climate and global carbon cycle (Rodenbeck et al 2018). El Niño and La Niña events occur every 2-7 years, when sea surface The GloSea5 seasonal forecasting system is based on a coupled climate model, which includes simulation of the land surface and vegetation in its component model JULES (Joint UK Land Environment Simulator, Best et al 2011, Clark et al 2011. Although GloSea5 does not include a coupled carbon cycle, NPP is calculated internally by JULES. Fire risk indices can also be calculated offline, based on GloSea5 forecast model output. We are therefore able to examine how well the same seasonal forecasting system used by Betts et al ( , 2018 can directly model separate elements of the carbon cycle. In this paper we assess the skill of GloSea5 for forecasting NPP and the McArthur Forest Fire Danger Index (McArthur 1966, Luke andMcArthur 1978) in the Tropics, and explain this with reference to the skill in the underlying meteorological variables. We also investigate a simple fire index based only on temperature and relative humidity, to understand the impact of using fewer and simpler quantities in the calculation of fire risk on seasonal time scales.
We describe the data sets, methods and calculations we use in our analysis in section 2. Section 3 describes our results, firstly validating our observation-based data, then considering seasonal forecast skill. We discuss our conclusions in section 4.

Data and analysis methods
In this section we describe how the fire risk indices and NPP are defined, and how they are calculated for our observation-based reference values. We then describe the GloSea5 seasonal forecasting system, and how we use that to calculate hindcast-based values of the fire indices and NPP.

Fire risk indices
Fire indices are used operationally in many fire-prone regions to help risk reduction and planning. A wide variety of indices have been developed to quantify fire risk, recently reviewed by de Groot et al (2015). Although they are usually tuned to specific regions, some have been used successfully in global-scale climate simulations (e.g. Betts et al 2015, Burton et al 2018), including the two we consider here: the McArthur index and the Angström index.
The McArthur Forest Fire Danger Index Mark 5 was first developed for Australia, and is calculated on a daily basis, integrating information on past rainfall leading up to each day (Noble et al 1980 Keetch and Byram (1968) drought index f drt has a maximum limit of 10 (Sirakoff 1985), and is given by where N is the number of days since the last day with rainfall, and R is the amount of rainfall on that day in mm. The restore amount, a restore , is the amount of water needed to restore the water content of the soil to field capacity, given in the range 0-200 mm. Given the difficulty in obtaining global soil moisture observations, we remove the soil moisture dependence entirely, by setting a restore to a constant value of 120 mm, as used in Golding and Betts (2008). We have tested the effect of this simplification on correlations of seasonal means in appendix A.
The McArthur fire index is used particularly in Australia, South Africa and Spain (de Groot et al 2015), but it has also been shown to strongly correlate with satellite estimates of actual fire occurrence in the Amazon (Hoffmann et al 2003), and Golding and Betts (2008) used it in their study of Amazonian fire risk under future climates. Note that the McArthur index is defined such that higher values indicate greater risk.
The Angström index has a much simpler definition: where RH % is daily mean relative humidity as a percentage, and T C is the daily mean temperature in°C. , following Eastaugh et al (2012); other constants are also in use, e.g. Skvarenina et al (2004) and Holsten et al (2013). Since this equation is linear, it can also be used with monthly or seasonal mean relative humidity and temperature values. Note that smaller values of the Angström index indicate greater fire risk.
The Angström index was developed in the first half of the Twentieth Century in Sweden, and has been used throughout the Scandinavian Peninsula (Hamadeh et al 2017), but has had success in other countries such as Slovakia (Skvarenina et al 2004), Taiwan (Lin 1995), Germany (Holsten et al 2013), Austria (Arpaci et al 2013) and Brazil (Alves White et al 2013).
We calculate observation-based reference values of both fire indices, using the WFDEI reanalysis data set (WATCH Forcing Data methodology applied to ERA-Interim, Weedon et al (2011, 2014. This covers land points only, on a 0.5°grid. We calculate the reference Angström index from monthly mean temperature and relative humidity. We calculate the reference McArthur index using daily maximum temperature, daily minimum relative humidity, daily mean windspeed, and daily precipitation. We use a near-zero threshold for daily precipitation (10 −8 mm), but we have also tested a higher level (5 mm) following Keetch and Byram (1968). The seasonal-mean results from using either threshold using the WFDEI data are very highly correlated.
Although these fire indices do not account for available fuel in their quantification of fire risk, we mask out gridcells where the fraction of bare soil is 0.5, following Gilham (2014). This removes regions without substantial amounts of vegetation to burn.
It is important to understand how the fire risk indices we consider here relate to actual fire occurrence. In principle, there could be significant differences: the fire indices do not explicitly account for fuel availability or chance of ignition, for example. These indices are based solely on meteorological factors, as opposed to human factors. We use observations of burnt area to quantify fire occurrence, using v4.1s of the Global Fire Emissions Database (GFED, Giglio et al 2013, which includes an experimental adjustment for small fires (Randerson et al 2012. The GFED data starts in 1997, later than the WFDEI and seasonal forecast data sets we use.

Net primary productivity
Net primary productivity is calculated from the carbon assimilated by photosynthesis minus the carbon lost by respiration. This quantity is not directly observable. At the site level, eddy correlation techniques can be used to measure net ecosystem exchange, which incorporates both NPP and soil respiration. Field-based measurements can give accurate NPP data on a small scale, but it is difficult to estimate NPP on larger scales due to sparse observation networks (Pachavo and Murwira 2014). Satellite products exist (e.g. MODIS NPP, Zhao et al 2006, Heinsch et al 2006, but these are themselves based on models, and they can be contaminated by cloud coverparticularly relevant in the Tropics. The MODIS NPP data is only available from 2000, which would give too short a period to robustly assess the skill of seasonal forecasts. To avoid these problems, and for greater internal consistency, we calculate our observation-based reference values of NPP using JULES driven by the WFDEI reanalysis data. In JULES, NPP depends strongly on soil moisture (Clark et al 2011, Harper et al in prep). Soil moisture depends on evaporation and transpiration, which in turn depend on vegetation cover and atmospheric conditions (in particular, temperature and humidity, with wind also being a factor). NPP in JULES also has a strong dependence on temperature (Clark et al 2011), via both photosynthesis and respiration.
Modelled primary productivity in JULES has been evaluated against a number of observational datasets, and is shown to successfully reproduce interannual variability of GPP on global scales, while on a regional scale GPP in the tropics was biased to higher values than observed (Slevin et al 2017).
The JULES configuration used for our WFDEI-driven simulations is based on 'JULES-C' as used in the Global Carbon Budget annual assessments (le Quéré) and assessed for its ENSO response characteristics (Bastos et al 2018). In the present application the dynamic vegetation model was switched off and the vegetation cover prescribed using the International Geosphere and Biosphere Programme (IGBP) land classification dataset, which was processed as part of the WFDEI project. The soil ancillary was generated by the Central Ancillary Program (Dharssi et al 2009) on the WFDEI grid, using the Brooks and Corey (1964) parameterisation. The topographic index data was from Hydro1k (EROS Archive 2017). The model was forced with global CO 2 from NOAA. 4

Seasonal forecasts
We use data from the GloSea5 seasonal forecast system (MacLachlan 2015) in its Global different contexts and applications, of particular relevance here is that GloSea5 has very high levels of skill in forecasting ENSO (MacLachlan 2015), and in forecasting rainfall in the tropics (Scaife et al , 2019. Note that the version of JULES used within GloSea5 is less recent than the one we use to calculate observation-based NPP (see previous subsection).
Calculating the Angström fire index from the hindcast data is reasonably straightforward. As already discussed, the Angström index is linear, so its seasonal mean values can be calculated directly from monthly mean forecast output. In contrast, the McArthur index can only be approximated by the forecast data. Although daily maximum temperature and daily mean wind speed were available, the daily minimum relative humidity was not stored when the forecasts were made, so we have used the daily mean in that case. A second difficulty is in the calculation of the drought factor: the number of days since there has been any rainfall could be longer than the forecast lead time in some locations (e.g. before 25th Oct, 1st Nov or 9th Nov for the DJF forecasts). Those cases are masked out of our calculations.
To make an initial assessment of the forecast skill, we will be using the simple interannual Pearson correlation between the seasonal mean, ensemble mean hindcast data, and the reference data derived from observations. Using the ensemble mean ensures that we are maximising the predictable signal from the ensemble; having skill in the ensemble mean is a necessary first step before examining more complex skill metrics based on the distribution of ensemble members.
We are not applying any bias correction to any of the hindcast data. While the Pearson correlation is not sensitive to overall seasonal-mean biases, the fact that there are many absolute thresholds in the fire index definitions means that model biases could still affect the results. We are also not detrending any of the data when assessing its skill. While it is often informative to understand how much of a skillful forecast comes from a trend versus other processes, in our case it is just as relevant to know that the model is able to reproduce an observed trend. Studies detailing the ability of a seasonal forecast to produce a useful climate service for a particular user would require more detailed consideration of biases and trends, as well as probabilistic skill assessment. For this initial study, of how well the initialised climate model on its own can capture the processes necessary for forecasting relevant components of the carbon cycle, we simply test the correlation.

Results
We first show how tropical NPP and fire risk behave under El Niño and La Niña conditions in our observationbased data sets. This provides important context for the subsequent subsections on their performance within the GloSea5 hindcasts. For consistency, all plots cover the same 20-year period and two seasons as the hindcast (1992/93-2011/12, DJF and JJA), unless noted otherwise. Figure 1(a) shows, for meteorological context, the mean anomalies of seasonal rainfall totals under El Niño and La Niña conditions, based on the Oceanic Niño Index 5 (ONI): A season is labelled as being El Niño if its seasonal mean value is >0.5 K; and as La Niña if its ONI <−0.5 K. We classify DJFs and JJAs separately according to their own ONI values. There are clear regions of enhanced or reduced precipitation across the Tropics, with impacts in different regions at different times of the year, reflecting the changes in atmospheric circulation induced by the ENSO events.
3.1. NPP and fire index relationship with El Niño figure 1(b) shows the response of NPP to El Niño and La Niña events in DJF and JJA. There are regions of both reduced and enhanced productivity in each case. The response matches that of precipitation in water-limited regions of northern South America in DJF, and in Africa in both seasons. The NPP signal in northeastern Brazil in JJA is not reflected in the precipitation response, but is reflected in variables that affect soil moisture via evaporation. There is no NPP signal in Indonesia in JJA despite the clear precipitation response, likely due to vegetation in the region not being water limited.
Similar composite maps are shown for the McArthur and Angström fire indices in figure 1(c) and (d). These follow the precipitation response more closely, with enhanced fire risk in northern South America, southern Africa and northwestern Australia in DJF El Niños, and in Indonesia and eastern India in JJA El Niños. As with NPP, there are some additional responses not seen in the precipitation, such as in the Sahel, the Horn of Africa and India in DJF, and northern Australia in JJA. figure 2(a) shows the ENSO response in the observed burnt area data. For consistency with our other results, we only use the period up to the end of the GloSea5 hindcast data set, i.e. the 15 years of 1997/98 to 2011/12 for DJF, and 1997-2011 for JJA. As the GFED data set is based on very high resolution observations of individual fires, this data is much more spatially varied than the meteorology-based data shown previously (all our figures use data regridded to match the GloSea5 data). Nevertheless, there are clear similarities with the fire risk indices, particularly in northern South America, the Sahel and northwest Australia in DJF. It is important to note that the two seasons we focus on here do not necessarily correspond to the peak fire seasons around the world; for completeness, we show the burnt area ENSO response in the other two seasons, and the climatologies of all four seasons, in appendix B.

Fire index relationship with burnt area
To clarify the relationship between burnt area and the fire risk indices, figure 2(b) maps the correlation of their seasonal mean time series. This shows that, in almost all the regions where the fire indices respond strongly to ENSO, they are well correlated with observed burnt area (although a noticeable exception is northern Australia in JJA). This gives us confidence that the fire indices give a reasonable estimate of fire risk.
One of the most important features seen in figure 1(c), (d) and figure 2 is that there is very little difference between the results for the McArthur and Angström fire indices: they both perform equally well/poorly in terms of following the interannual variability of burnt areas in the tropics. This is backed up by plotting the correlation between the two fire indices themselves (figure 3): they are very highly (anti-)correlated, with values <r 0.8 in most places for both seasons, over the 20-year period we consider here. The additional complexity in modelling fire risk that goes into the McArthur index does not seem important for most regions of the tropics for forecasting seasonal means.

Seasonal forecasting skill
In the following figures, we show seasonal forecast skill maps, for NPP (Figure 4) and the McArthur and Angström fire indices (figure 5). Figure 6 shows the skill of their underpinning meteorological variables: temperature, precipitation, relative humidity and wind speed.

NPP and its components
In DJF, there is skill in NPP forecasts (figure 4) in equatorial South America, Africa, India and the Indochinese Peninsula. However, air temperature, precipitation, relative humidity and wind speed all show greater skill over much larger areas of South America in DJF (figure 6), particularly in the region from the Venezuelan coast south to the Amazon river. This is a good demonstration of how skill in the underlying variables does not map directly on to skill in the compound variables. In India and Indochina, the high correlation in NPP seems more related to temperature, with Indochina in particular not showing widespread skill in relative humidity or precipitation. These are also regions where NPP does not show a clear ENSO response in DJF.
In JJA, the correlation between reanalysis-forced NPP and GloSea NPP is high in northern Venezuela, northeastern Brazil, eastern Africa and northern Australia. These regions show a clear response to El Niño in JJA ( figure 1(b)). The NPP skill shown could have useful applications, such as for agricultural production in northeast Brazil. In this region, the main contributor to the high skill is likely to be the depletion in soil moisture due to evapotranspiration, affected by relative humidity and wind speed, which both show skill in this area.

Fire indices and their components
There are large areas of skill in both DJF and JJA for the McArthur fire index, with very similar patterns seen for the Angström index skill ( figure 5). The areas of skill largely correspond to where the fire indices in these seasons  This covers the same 20 year periods as the GloSea5 hindcast; a correlation must be > r 0.44 | | to be significantly different to zero according to a Fisher z-test at the 5% level, and this test is passed almost everywhere. respond strongly to ENSO (figure 1(c) and d). Comparing with the skill of the meteorological variables (figure 6), it is clear that the skill patterns are very similar to those of relative humidity and, to a lesser extent, wind speed and precipitation. The temperature skill patterns differ, but temperature is more skillful generally, over broader areas, partly due to its trend. This points to the key processes needed for forecasting fire risk being simply a measure of dry and hot conditions. Finally, we estimate the skill of GloSea5 in forecasting burnt area, using a fire index or a single meteorological variable as a linear predictor (i.e. simply the correlation between a variable from the GloSea5 hindcast and the observation-based burnt area data). Figure 7 shows the Pearson correlation for each of these variables with burnt area. Here, a strong positive or negative correlation indicates that the variable could be used as a skillful predictor.
Although the skill is clearly reduced compared to forecasting the fire indices or meteorological variables themselves, there are nevertheless some coherent large-scale regions that suggest useful levels of skill: northern South America, southern Indochinese Peninsula, northwestern Australia in DJF; northeastern Brazil, east Africa and Indonesia in JJA. It also appears that, to a large extent, just using a single meteorological variable could provide skillful forecasts of fire risk in some regions.
As expected from our examination of the observation-based data, the additional complexity of the McArthur index does not significantly improve the forecast skill, either of the fire indices themselves, or of actual burnt area. . Correlation between ensemble-mean data from the GloSea5 hindcast, and observation-based data for NPP, for DJF (left) and JJA (right). A green contour marks = r 0.44; | | correlations larger than this are significantly different to zero at the 5% level, based on a Fisher z-test.

Discussion and conclusions
This study has examined the seasonal forecasting skill of two key processes that contribute to the interannual variability of global CO 2 levels: fire risk and net primary productivity. While we expect that hybrid dynamicalstatistical forecasts (e.g. based on ENSO, similar to Betts et al , 2018 will still be the best approach to forecasting global CO 2 , our investigation is useful for highlighting underlying capabilities in the climate model used for our seasonal forecasts. We highlight where the forecasts do and do not have skill, in different regions and at different times of the year, in order to determine whether the skill is in the key regions and seasons that contribute to the CO 2 variability. The regional distribution of skill can also inform approaches to seasonal forecasting, and, in principle, the development of climate services for agriculture and fire risk management.
The seasonal forecast skill for the NPP modelled explicitly in GloSea5 is significant for large areas of the Tropics, including where NPP responds strongly to ENSO. Some of the regions where it is skillful are well-placed for forecasting crop yield, such as northeast Brazil.
The fire indices we have considered also show significant skill across much of the Tropics. However, we have also demonstrated, in observations as well as in terms of seasonal forecast skill, that on seasonal-mean time scales there is no benefit to using the detailed McArthur fire index over a simple indicator of hot and dry conditionsthe Angström fire index. Any benefits from having a more detailed, carefully-calibrated fire index, as might be seen on weather-forecasting time scales, are lost when averaging over a season. The resulting skill picks out largescale features in common across both indices: areas where it becomes anomalously dry and hot at similar times. We have demonstrated that there are several regions where burnt area itself could be forecast skilfully, based either on a fire risk index, or perhaps even on a single meteorological variable.
A third process, which is not considered here, but is also an important contributor to the interannual variability of terrestrial carbon emissions, is heterotrophic respiration. This is currently modelled within the GloSea5 system but not available as an output variable. Model heterotrophic respiration includes a factor that depends explicitly on soil temperature, and a factor that depends explicitly on soil moisture (Clark et al 2011). An interesting extension to this work would be to look at whether the GloSea5 system can skilfully predict these two factors, and how much of this skill can be captured for seasonal means using a simple combination of meteorological variables such as air temperature, precipitation and wind speed.
Here, we have focused on the terrestrial components of the carbon cycle that are known to dominate the interannual variability of atmospheric CO 2 (Jones et al 2001). However, the ocean plays a smaller, but still Figure 7. Correlation between the observed burnt fractions using the ensemble-mean fire indices and meteorological variables (as labelled) from the GloSea5 hindcast. The correlation for the Angström index has had its sign switched, for consistency with the McArthur index. Significance contours are not drawn, but as this only uses 15 years of data, a nominal threshold for significance at the 5% level would be > r 0.51 | | . significant, opposing response to the land for a typical El Niño event, of approximately 25% of the land magnitude. A logical next step towards developing a physically-based seasonal atmospheric forecast would be to assess the skill of both terrestrial heterotrophic respiration (as discussed above) and the air-sea flux of carbon. Ocean biogeochemistry is not currently included in the GloSea5 forecast system, although the related UK Earth System Model (UKESM, Sellar et al 2019) does include an ocean biogeochemistry model, MEDUSA. However, it is computationally very expensive, and including it in forecasts may not add predictive skill. This again implies that a more parsimonious approach, through understanding the seasonal drivers of ocean carbon exchange, is likely to yield better results. Our results are a clear demonstration of the different kinds of analysis required when considering seasonal climate prediction, compared to weather forecasting studies, which require confidence at forecasting individual fire events; or climate projection timescales, where longer-term feedbacks are crucial. When forecasting seasonal means, simplicity is important. An impacts model that requires driving data from many input variables, at high temporal resolution, is unlikely to result in any benefit in terms of skill; much simpler metrics based on one or two variables on monthly-to-seasonal timescales are likely to be a more efficient way of achieving skillful forecasts. Palin et al (2016) and Bett et al (2019) have demonstrated this point in the very different contexts of seasonal forecasting for transport and energy sector impacts in Europe.
This also applies to interactive models, like JULES' fire model INFERNO ( Development of climate services, for crop yield or fire risk, would require more detailed assessments, considering bias correction and probabilistic forecast calibration, similar to the study by Bedia et al (2018) on seasonal forecasting fire risk in the Mediterranean. Although they used the Canadian Fire Weather Index, which is similar in complexity to the McArthur index used here, they also found that their regions of good skill were closely linked to skill in forecasting relative humidity. Being able to use simpler impacts metrics could greatly simplify the development and production of seasonal climate services for forecasting fire risk. However, it could be the case that, operationally, the most important quantity to forecast is not the seasonal mean, but rather the number of high fire risk days within the season, or some other threshold-based quantity. This should be determined in conjunction with the relevant stakeholder as part of the co-development of the forecast service, as demonstrated by Turco et al (2019): they describe the development of a climate service prototype for probabilistic forecasts of the burnt area in Catalonia being above a user-defined threshold.
Our study is intended to be a starting point for investigating model capability for seasonal forecasting of carbon cycle components. There are of course many limitations: The NPP reference data are reanalysis-driven model runs; it was necessary to approximate the fire index calculations due to limitations in available data; and the number of years in the hindcast results in large uncertainties in the correlation skill. However, we have demonstrated that GloSea5 does have significant skill in forecasting the year-to-year variation in NPP and fire risk, which are important components in the carbon cycle processes that contribute towards global CO 2 variability. tested using an explicit soil moisture (θ) dependence by calculating a restore as the amount of water needed to restore the top 1 m of soil to field capacity (q FC , at −0.01 MPa), as a fraction of the difference between field capacity and soil moisture wilting point (q W , at −1.5 MPa), rescaled to be in the interval 0-200 mm: We used the daily soil moisture from the same WFDEI-driven JULES simulations we ran for our reference NPP calculations, to calculate a new McArthur fire index data set for testing.
Despite the restore amount showing significant regional variation ( Figure A1), the seasonal mean McArthur index based on varying soil moisture, and the simplified index using a constant restore amount, are very highly correlated (figure A2).
Appendix B. Burnt area seasonal climatologies throughout the year Figure B1 shows the climatological seasonal mean fields for burnt fraction of gridcells, covering all four seasons: DJF and JJA as in the body of this paper, plus March-April-May (MAM) and September-October-November (SON). This illustrates the regional distribution of the peak fire season throughout the year. Figure B2 shows, for completeness, the response of the burnt fractions to ENSO, as in figure 2(a), but for MAM and SON.  ORCID iDs