Estimating grassland curing with remotely sensed data

Wildfire can become a catastrophic natural hazard, especially during dry summer seasons in Australia. Severity is influenced by various meteorological, geographical, and fuel characteristics. Modified Mark 4 McArthur’s Grassland Fire Danger Index (GFDI) is a commonly used approach to determine the fire danger level in grassland ecosystems. The degree of curing (DOC, i.e. proportion of dead material) of the grass is one key ingredient in determining the fire danger. It is difficult to collect accurate DOC information in the field, and therefore ground-observed measurements are rather limited. In this study, we explore the possibility of whether adding satellite-observed data responding to vegetation water content (vegetation optical depth, VOD) will improve DOC prediction when compared with the existing satellite-observed data responding to DOC prediction models based on vegetation greenness (normalised difference vegetation index, NDVI). First, statistically significant relationships are established between selected ground-observed DOC and satellite-observed vegetation datasets (NDVI and VOD) with an r2 up to 0.67. DOC levels estimated using satellite observations were then evaluated using field measurements with an r2 of 0.44 to 0.55. Results suggest that VOD-based DOC estimation can reasonably reproduce ground-based observations in space and time and is comparable to the existing NDVI-based DOC estimation models.


Introduction
Wildfire can be responsible for major environmental damage or changes to ecosystems (Cobb et al., 2016;Gazzard et al., 25 2016;Mistry et al., 2016). One of the important components in determining the severity of wildfire is fuel availability. Wildland fuels can vary considerably, both spatially and temporally (Stambaugh et al., 2011). Various interpretations and characterisations of fuel have been made in past studies as a key contribution to assessing wildfire potential (Hudec and Peterson, 2012;Jurdao et al., 2012;Sharples et al., 2009b;Stambaugh et al., 2011;Yebra et al., 2013). Fuel can also be quantified by its age or time since last fire (Bradstock et al., 2010). 30 In this study, we focus on the availability of combustible fuel in the aboveground biomass in grassland ecosystems; this fuel availability metric is referred to as the degree of curing (DOC). The DOC is the percentage of dead material in a grassland fuel bed; 100 % indicates a fully cured (dead) grassland fuel complex. The DOC has a direct influence on wildfire development in grasslands, hence, it is an important input for fire danger indices and fire spread models, such as the McArthur Grassland Fire Danger Index (GFDI) (Gill et al., 2010) and the CSIRO grassland fire spread model Kidnie et al., 2015). 35 Generally, fires are unable to spread across grasslands that are less than 50 % cured . Though this lower limit has been revised, since a more recent study demonstrated that fire can spread in grassland with DOC as low as 20 % . In climatological studies, DOC is often assumed to be 100 % (Pitman et al., 2007). This leads to an overestimation of areas experiencing high levels of fire danger and hence provides only a weak indication of where to focus resources of fire agencies. An accurate, spatially and temporally explicit, estimate of DOC would provide more useful guidance to these agencies.
Measuring DOC in the field is a tedious and expensive task, especially when an accurate assessment of curing is required. Anderson et al. (2011) suggested that current methods for measuring DOC still present problems. The visual assessment 5 method, which relies on field observers to estimate the general curing value based on their expertise and a visual guide, is subjective and can often be unrepresentative of DOC of the entire area. Destructive sampling approaches can provide accurate field based observation, but is a labour intensive task. Thus, Anderson et al. (2011) offered a simple field based method utilising a levy rod, based on the modified point quadrant method of pasture assessment; the approach involves counting the number of live and dead touches on a thin steel rod that was driven into the ground. It was suggested that this approach can be 10 applied across Australia with higher accuracy than current visual assessment methods .
Apart from ground measurement, DOC can be estimated using satellite remotely-sensed data, but it is limited by the satellite sensors' capability, e.g. spatial resolution and various atmospheric interferences. Dilley et al. (2004) established a relationship between curing and Normalised Difference Vegetation Index (NDVI, a proxy of vegetation canopy greenness) by estimating live fuel moisture content from NDVI and relating it to curing via an exponential function using a finite difference  Marquardt method (Dilley et al., 2004;Rouse et al., 1973). Newnham et al. (2011) showed that estimation of curing using a relative greenness (RG) approach that was based on NDVI distribution provided more accurate estimation of curing than a direct linear regression between curing and NDVI. (Chladil and Nunez, 1995) used curing derived from a soil dryness index model and NDVI to predict soil and fuel moisture content. Various optical based vegetation indices computed from remote sensing reflectance products can also be developed into a satellite based model integrated with ground observations to predict 20 curing Turner et al., 2011). These methods, though vastly different in their approaches, achieved good results for their set objectives, but they tend to focus on particular applications. It should also be noted that optical based remote sensing products, including NDVI, are affected by cloud cover and aerosols. Some studies explicitly acknowledge challenges presented by cloud effects and when there are both forest and water bodies in the same NDVI pixel, which results in an erroneous grassland interpretation (Allan et al., 2003;Chladil and Nunez, 1995). However, if appropriately detailed aerosol 25 data is available, atmospheric correction can mitigate the aerosol effect on NDVI.
Currently, there are satellite based DOC products over Australia provided by Bureau of Meteorology. The products have 500 m spatial resolution and 8 day temporal resolution, and are based on two past studies Newnham et al., 2010). There are five separate satellite based DOC models; four are from Newnham et al. (2010) and one is from Martin et al. (2015). All satellite based DOC models here are based on optical and near-infrared wavelength bands. We would like to 30 investigate whether including a recent passive microwave based satellite product can improve the DOC estimation over Australia or not.
A passive microwave based remote sensing vegetation product, referred to as Vegetation Optical Depth (VOD), has been developed recently (Meesters et al., 2005). VOD is primarily sensitive to vegetation water content, including both leafy and woody components (Guglielmetti et al., 2007;Jackson and Schmugge, 1991;Kerr and Njoku, 1990). Unlike the traditional 35 optical based vegetation indices, such as NDVI, VOD is minimally influenced by the atmospheric conditions due to its longer wavelength and stronger penetration capacity (Jones et al., 2009). However, it has a coarser spatial resolution (0.1°) in comparison with optical based products, which is a consequence of the low energy microwave emissions from the Earth's surface. It has been demonstrated that VOD can capture the changes in vegetation water content over different land cover types at the global scale, including grassland, cropland, savannas, tropical forests, and boreal forests (Liu et al., 2013a(Liu et al., , 2013b(Liu et al., , 2015. 40 A recent study on estimating fire severity based on VOD changes also demonstrated that VOD has high correlation with aboveground biomass in Australian tropical savannahs that includes large area of northern grassland with an r 2 of 0.96 (Chen et al., 2018). Also, NDVI and VOD provide complementary information and can comprehensively characterise vegetation dynamics when combined (Andela et al., 2013).
There are two objectives of this study. The first is to explore the possibility of whether adding the VOD (responding to vegetation water content) will improve DOC prediction when compared with existing NDVI (responding to vegetation greenness) based DOC prediction model. The second is to implement the satellite based DOC estimation (both our and existing 5 models) into GFDI to investigate whether better fire severity predictions can be achieved in grassland environments.

Satellite based Products
The NDVI dataset used here is derived from the Moderate Resolution Imaging Spectroradiometer (MODIS) 8 day surface reflectance product (MOD09A1) on-board the Terra satellite (Vermote and Vermeulen, 1999). The MOD09A1 product used 10 here for computing NDVI is an 8 day product, which has less noise than the daily product, and its spatial resolution is 0.005° (~500 m). This product is obtained from Remote Sensing at National Computational Infrastructure (NCI) MODIS Land Product for Australia website (Paget and King, 2008). It is produced from original tiles provide by the United States Geological Survey (USGS) for Australia with a starting date from February 18th, 2000.
The 8 day NDVI data is derived from the MODIS reflectance dataset using the following Eq. (1): 15 where ρ 1 and ρ 2 are spectral reflectance measurements obtained from the visible (red) and near-infrared regions, respectively.
During the conversion, to ensure the quality of data, only pixels with ideal quality in all bands and a view angle zenith of less than 60° were kept for the analysis, as suggested by Newnham et al. (2011). The spatial resolution is kept at 0.005°.
The VOD dataset used here is retrieved from the Advanced Microwave Scanning Radiometer -Earth Observing System 20 (AMSR-E) and derived using the Land Parameter Retrieval Model (LPRM) approach from which soil moisture and VOD are retrieved simultaneously (Meesters et al., 2005;Owe et al., 2001). Several assumptions are made in the LPRM approach, including: canopy surface temperature equal to soil surface temperature, a constant single scattering albedo, same vegetation parameters for both Horizontal and Vertical polarizations, and minimal effect of surface roughness (Meesters et al., 2005;Owe et al., 2001). Uncertainties in soil moisture and VOD retrievals are expected with these assumptions. The evaluation of LPRM 25 soil moisture over Australia showed that the temporal patterns of satellite-based and in situ soil moisture agree very well (Draper et al., 2009;Gevaert et al., 2016). This agreement suggests a reasonable separation of temporal patterns of soil moisture and VOD, while uncertainties may exist in the absolute magnitudes of these two variables.
VOD has a spatial resolution of 0.1° (~10 km) and nearly daily temporal resolution (Parinussa et al., 2014). The time period covered by AMSR-E is from 2 June 2002 to 3 October 2011, but we used a 9 year range from 4 July 2002 to 26 June 2011 in 30 our analysis by excluding the beginning and the end of the AMSR-E records. The same time period is also used for the MOD09A1 NDVI dataset. VOD data, which has a near daily temporal resolution, is converted to an 8 day average product to reduce noise and ensure complete coverage over Australia (10° S to 45° S and 110° E to 160° E) per temporal interval. The grid cells with radio frequency interference (RFI) are excluded from our analysis. RFI is caused by man-made transmitters, such as radars and 35 wireless communications. These transmitters can be operated in the same frequency range as passive microwave observation, including VOD. Thus, natural signals captured by passive microwave observations are sometimes contaminated with RFI (Nijs et al., 2015). Vegetation types that are present within the VOD (0.1°) and NDVI (0.005°) pixel influence the differences in VOD and NDVI behaviour. The difference between VOD and NDVI spatial resolution can be clearly seen in the example 2° by 2° spatial maps around the Silent Grove area (Fig. 1). We use both VOD and NDVI together to combine their strengths.
The global 0.05° land cover map based on the MODIS MCD12C1 product is used for classifying the dominant land cover type within each VOD pixel. The land cover classification system is as proposed by the University of Maryland (UMD scheme) (Hansen et al., 2000). Figure 2 shows the 0.05° land cover type map of Australia, with observed curing sites marked with 5 crosses. The land cover map used here is from year 2010.
A burned area product from MODIS is acquired for further evaluation of the recalculated GFDI from satellite based DOC results. The monthly archived MODIS burned area map reprojected for Australia is obtained from Remote Sensing at the NCI site (Paget and King, 2008). There are two separate MODIS burned area products: the MCD45A1 and the MCD64A1. The MCD64A1 burned area product is preferred over MCD45A1, since it was proven to be more accurate (Andela and van der 10 Werf, 2014;Padilla et al., 2015;Ruiz et al., 2014). Its spatial specification is exactly the same as the MODIS reflectance dataset, with temporal availability from August 2000 onwards. To ensure high quality of the burned pixels, only pixels with the valid data flag from the provided quality control file are included in the analysis. Over 99 % of pixels from mid-2002 to mid-2011 are classified as unburned. To reduce the number of prescribed burned and other low power anomalies detected by the burned area product, a fire radiative power (FRP, unit: MW) from MODIS active fire product (MCD14ML) is used to 15 mask out low severity fires.

Ground Based Curing Observations
The observed grassland DOC data was provided by Bushfire and Natural Hazards Cooperative Research Centre and its partner agencies [Project reference: http://www.bushfirecrc.com/projects/a14/grassland-curing]. The observed data were collected from several sites across Australia and New Zealand, ranging from August 2005 to March 2009, usually during summer (Fig.  20 2). Selected sites were intended to represent broad coverage of major grassland types. Note that the number of locations and samples taken were highly dependent on the availability of field observers from fire management agencies; data were collected with inconsistent interval between data collection dates . Three types of data collection approaches were used: visual estimation, levy rod method, and destructive sampling. Due to the number and availability of data as well as their accuracy, only observed DOC from the levy rod method was used in this study. Anderson et al. (2011) and Newnham et al. 25 (2011) state that the levy rod measurement is reliable with RMSE of 13.5 % and a bias of less than 1 % when compared with destructive sampling.
To identify robust relationships between the site observed and remotely sensed DOC, a number of site criteria must be met.
Sites meeting these criteria were used for calibration of the VOD and NDVI satellite based DOC models, while all of the valid records were used for evaluation. One major factor in deciding the site selection is the land cover properties of the observed 30 DOC site. The 0.05° land cover type map (MCD12C1) is used for classifying the site location land cover (Hansen et al. 2000).
Since 0.1° VOD pixel is covered by 2 by 2 0.05° land cover pixels, the corresponding 2 by 2 pixels of land cover type for each observed curing site can be acquired. The land cover type and homogeneity of each observed DOC site can then be determined, where the site is considered to have a homogeneous land cover only if all four land cover pixels corresponding to the VOD pixel are the same. In case of a site with heterogeneous land cover type, the dominant land cover with the most pixels out of 35 four will be considered as the representative land cover. All observed DOC sites can be categorised into the following land cover types: evergreen broadleaf forest, open shrubland, savannas, woody savannas, grasslands, croplands, and urban.
According to the land cover information for each 0.1° VOD pixel, sites identified as evergreen broadleaf forest pixels were removed from the analysis. There are three out of 37 sites situated in the evergreen broadleaf forest, which are Darnum, VIC (mixed grass), Simcocks, WA (improved pasture), and Neerim South, VIC (mixed grass). Even though actual locations of all 40 observed curing sites were in grassy areas, the VOD signal is a mixture of grassland and forest when the sites are surrounded by dense forests within the same 0.1° pixel.
All sites were also examined to ensure a negative correlation between VOD and the in situ DOC data. That is, since VOD is a proxy for water content in above ground biomass, an overall negative correlation between VOD and curing is expected. If this is not the case, then there is likely some other activity within the 0.1° pixel that disrupts this basic relationship; this effect was found in three sites (Durran Durra, NSW (native grass), Monaro, ACT (improved pasture), and Parry Lagoons, WA (native grass)). The sites without the negative correlation with VOD also had no correlation with NDVI suggesting other land cover 5 types where dominating the signal. Thus, six out of 37 sites are excluded from the analysis.
In addition, there are eight sites (Umbigong, ACT (native grass), Kilcunda, VIC (improved pasture), Tooradin, VIC (improved pasture), Tooradin North, VIC (improved pasture), Caldermeade Park, VIC (improved pasture), Kaduna Park, VIC (improved pasture), Hobart Airport, TAS (native grass), and Jerona, QLD (native grass)) in which VOD data are not available. Most of these are due to sites being located near the coast or a large body of water, where the VOD signal is strongly influenced by the 10 water itself. With the remaining 23 out of 37 sites, several site selection criteria were applied for the calibration phase. The criterion used here to maintain consistency in observation time series requires sites to have at least eight consecutive records, where records are considered consecutive when they are separated by no more than 15 days. Only the consecutive series of records within the selected sites are included in the analysis for the calibration phase. This ensures that the derived model contains the temporal evolution of DOC within years. Only five out of 23 sites are retained for this group, containing a total 15 of 122 (out of 238 total) observations. The selected sites are: Majura, ACT (improved pasture), Tidbinbilla, ACT (mixed grass), Ballan, VIC (improved pasture), Murrayville 1, VIC (native grass), and Murrayville 2, VIC (improved pasture).
Multiple linear regression models of VOD and NDVI were then calibrated with the observed curing from the final selected sites.

Meteorological Datasets 20
To further assess the usability of the satellite based curing acquired from the VOD and NDVI model, the GFDI is computed.
Additional meteorological data needed for GFDI computation are dry bulb or maximum temperature, 3 pm relative humidity, maximum wind speed, and fuel load (Purton, 1982). Since fuel load is often set as a constant value of 0.45 kg m -2 (Sharples et al., 2009a), there are 3 remaining input datasets needed. These gridded, meteorological datasets are usually derived from the network of ground observation stations across Australia. The range for these datasets is from 4 July 2002 to 26 June 2011 to 25 exactly match with the VOD 9 year range. Both temperature and relative humidity datasets are acquired from the Australian Water Availability Project (AWAP) (Jones et al., 2009). Note that relative humidity is derived from vapour pressure and temperature data. These AWAP datasets have a 0.05° spatial and daily temporal resolution with a coverage region of 10° S to 44.5° S and 112° E to 156° E. For maximum wind speed data, the reanalysis maximum daily wind speed is computed from the ERA-Interim wind components dataset, acquired from the European Centre for Medium-Range Weather Forecasts (ECMWF) 30 (Dee et al., 2011). The reanalysis wind components dataset is available globally at approximately 0.8° spatial resolution at a 6 hour interval.

Developing VOD NDVI based dynamic DOC estimates
As indicated by past studies (Dilley et al., 2004;Peterson et al., 2008), NDVI has a significant relationship with live fuel 35 moisture content and DOC. In addition, the NDVI dataset has a spatial resolution of 0.005°. Thus, we investigate whether VOD adds any information to satellite based DOC beyond that embodied in the NDVI through the use of a multiple linear regression model. In addition, a past study suggested that a modified form of NDVI, referred to as relative greenness (RG), has stronger relationship with DOC than NDVI (Newnham et al., 2011). One of the two alternatives Newnham et al. (2011) proposed is the range based RG, which can be computed by the following Eq. (2) (2) where the NDVI min and NDVI max are the minimum and maximum NDVI value over a specified time range. Note that we did not attempt to compare the other purposed alternative, the spread based RG, but only range based (per pixel) RG. In Newnham et al. (2011), while range based RG performance is not as good as preferred spread based RG (r 2 = 0.62 and RMSE = 14.2 %), it is still better than plain NDVI (NDVI had r 2 = 0.50 and RMSE = 16.4 %, while 2.5 years range based has r 2 = 0.57 and 5 RMSE = 15.1 %). Note that while we cannot exactly reproduce 10 years time range Newnham et al. (2011) used, since our study time frame is 9 years, we tried various 2.5 years time ranges that overlapped with Newnham et al. (2011) study period.
However, given the same observed DOC data and NDVI dataset, we were unable to reproduce a result where the RG correlation with DOC is stronger than NDVI. Further analysis showed that the RG results were very sensitive to the selected time range for the computation, such that results were inconsistent with relatively small differences in the selected range. Due to this, RG 10 is not used in this study and NDVI is used directly in forming a multiple linear regression model to estimate satellite based curing. Some experimentation revealed that the VOD anomalies, computed from the difference between VOD and average VOD over a specified temporal range, yields the best correlation with DOC, but only if the VOD anomalies are computed from the range matching the in situ DOC observation range for each specific site. The range selection for computing VOD anomalies can be quite problematic, since it can heavily influence the correlation result, and no pattern could be found for 15 determining an appropriate VOD range for any other locations outside the observed DOC sites. Thus, we focus our analysis on using the absolute VOD value. The linear regression equation for curing and VOD correlation can be expressed as Eq. (3): where x 1 and x 2 are the intercept and slope of the relationship.
Utilising both VOD and NDVI datasets, the following multiple linear regression equation for estimating DOC can be expanded 20 from Eq. (3) as Eq. (4): where x 1 to x 4 are the intercept and coefficients of VOD, NDVI and the product of VOD and NDVI (interaction term), respectively. Using a stepwise regression, the calibrated final model with corresponding coefficients can be determined. The stepwise fit algorithm used here selects the significant terms with the lowest p-value, which is smaller than the entrance 25 tolerance, to be included in the model first. Next, the algorithm chooses the next most significant term that is still less than the entrance tolerance. This process is repeated until either there are no remaining significant terms or all terms are included in the final model (Draper and Smith, 1998). After the final model is calibrated, we evaluate the DOC model with all valid (23) observed DOC sites.

Comparing with existing DOC estimates 30
We acquire existing satellite based DOC products available from Bureau of Meteorology and compare their performance with our model. There are five models available, four are based on Newnham et al. (2010) and one is based on Martin et al. (2015) studies. We decided to test only one of Newnham's modelsthe one with the best overall RMSE (Method B), and Martin's model (MapVic). Both Method B and MapVic DOC models are as described in Eq. (5) where ρ 6 and ρ 7 are spectral reflectance band six and seven from MODIS reflectance dataset and GVMI is Global Vegetation Monitoring Index Newnham et al., 2010). GVMI can be calculated by:

Comparing GFDI dynamics using different DOC estimates
Several revisions of GFDI were made by past studies (Noble et al., 1980;Purton, 1982). The GFDI revision used in this paper 5 is modified Mark 4 GFDI, since it is the grassland fire danger assessing system that is generally being used by Bureau of Meteorology (Sharples et al., 2009b). Originally, the fire danger rating system was presented in a circular slide rule. A mathematical equation representation of modified Mark 4 GFDI was derived from the circular meter, and can be expressed as follows Eq. (8) (Purton, 1982): where Q is the fuel load (kg m -2 ), T max is the dry bulb or daily maximum temperature (° C), H 3pm is the daily relative humidity at 3 pm (%), V max is the daily maximum wind speed (km h -1 ), and f(C) is the curing factor. The curing factor can be calculated by Eq. (9): where C is the grassland DOC (%). 15 The GFDI is computed on the basis of meteorological input data and either a constant DOC at 100 % or satellite based dynamic curing values. These different GFDI datasets along with the burned area data (MCD64A1) can be used to examine the changes due to variable DOC spatially and temporally. By pairing up burned and unburned pixels with their associated GFDI pixel, we can assess the number of burned and unburned pixels for each GFDI severity level. Using histogram and receiver operating characteristic (ROC) analysis, the difference between original GFDI with constant DOC at 100 % and recalculated GFDI with 20 satellite based dynamic DOC can be assessed (DeLong et al. 1988;Zweig and Campbell 1993).

Comparing VOD-NDVI based and existing DOC estimates
Across all selected observed DOC sites (excluding the forest areas) from July 2002 to June 2011, the r 2 of VOD and NDVI is 0.52 with an RMSE of 0.11. Using the linear model, as described by Eq. (3), the DOC and VOD correlation result has a 25 significant relationship and an r 2 of 0.20 with RMSE of 20.80 %. The scatter plot showing the correlation between VOD, NDVI and combined VOD and NDVI terms are as shown in Fig. 3, while Fig. 4 shows the residual DOC unexplained by NDVI (differences between observed DOC and NDVI-based DOC) against VOD and combined VOD and NDVI terms.
However, at this level of r 2 , VOD alone is not reliable enough to estimate DOC, especially across Australia in general. The study of Newnham et al. (2011) indicated that NDVI alone can perform better at estimating DOC, with an r 2 of approximately 30 0.50 for a DOC and NDVI linear relationship. The combined explanatory power of NDVI and VOD is explored using a multiple linear regression analysis, as expressed by Eq. (4). The first final model includes the VOD and NDVI interaction ( 4 ) and NDVI ( 3 ) terms, and is as shown in Eq. (10). The calibrated r 2 for this model is 0.67 with RMSE of 13.40 %. VOD was excluded as a predictor in the first final model, as expressed in Eq. (10), because during the stepwise regression, when the NDVI and (VOD)(NDVI) terms are included as the first and second predictors, the VOD term does not contribute in improving 35 the final model prediction (i.e. p-value exceeds the acceptance threshold, preventing overfitting). When the NDVI term is excluded, the (VOD)(NDVI) term is included first, followed by the VOD term, as expressed in the second final model, shown in Eq. (11). The second final model has a calibrated r 2 of 0.54 and RMSE of 15.95 %. Table 1 shows the correlation results for both models.
The DOC models are then evaluated with all (23) valid observed DOC sites and independent (18) observed DOC sites (excludes 5 sites that were used in calibration). The evaluation results for the first model are also shown in evaluation and independent sites evaluation, respectively. While the evaluations resulted in degradation in model performance over the calibration in most cases, the independent evaluation of the second model has a slightly better evaluation performance.
These results can be compared to those obtained using existing remotely sensed DOC estimates which are also shown in Table   1. The MapVic DOC has a lower r 2 and higher RMSE, while the Method B DOC has a higher r 2 and lower RMSE when compared with both of our models in all sites evaluation. This indicates that Method B has the best evaluation among the three 10 models, while MapVic is the worst, and our models sit in the middle between the two. However, during an independent sites evaluation, both models with VOD have the worst performance (lowest r 2 and RMSE). This result is not entirely surprising as all the observations used here were also used in the calibration of Method B (Newnham et al., 2010), a subset is used in the calibration of our method, and MapVic was developed using an independent visual estimates dataset. That is, there is no independent data available for testing Method B, while both our method and MapVic are being tested against independent 15 data. While our first and second models do not have obvious advantages over one another, since the second model only performs better than the first in independent sites evaluation, we decided to pick the first model as our representative model for further comparison with the existing models from this point, since the terms in the first model were selected based on stepwise fit regression with none of our interference (we intentionally removed NDVI term in the second model before applying the stepwise fit regression). 20 Using the relationship between VOD, NDVI, and observed DOC from the first model, as stated in Eq. (10), we calculated satellite based DOC for Australia. Figure 5 presents maps of satellite based DOC data averaged over the summer periods To determine the amount of spatial variation in DOC across Australia, we computed the standard deviation of all valid DOC estimates across the continent within a single time step. All areas that are indicated as forests by the land cover type map are excluded from the analysis. The spatial variation time series can then be plotted for the available time period of mid-2002 to 30 mid-2011, as shown in Fig. 6. Note that the continental mean spatial DOC standard deviation is 20.39 %. This indicates that there is significant spatial variability in DOC that persists across all years, and contains a small seasonal component. For a normally distributed variable, 95 % of values would lie within two standard deviations, which is ±40.78 % in this case. Further analysis on DOC spatial standard deviation are as shown in Table 2. This includes seasonal, monthly, and land cover type spatial standard deviation of DOC. From both seasonal and monthly spatial standard deviation of DOC, it is shown that DOC 35 has the highest spatial variation during winter, which is especially true for northern Australia .
In addition, based on time series of satellite based curing data, Fig. 7 reveals the spatial distribution of standard deviations calculated for each pixel. It shows that most of the strong temporal variation occurs in the south, especially in the southeast and southwest of Australia. Several areas in the midcontinent that have unexpectedly high variation are likely due to rare inundation events. The continental mean temporal standard variation is at 11.88 %. Together, Fig. 6 and Fig. 7 show the 40 variability in DOC that will impact calculations of fire danger indices.

Comparing GFDI dynamics using different DOC estimates
The spatial plot for maximum summer recalculated GFDI from the DOC multiple linear regression model is shown in Fig. 8 The fire locations are marked with a red crosshair. White pixels are forest areas that were masked using the land cover map. 5 Overall, summer 2003 has 4.51 % more areas indicated as severe or higher GFDI than summer 2010. MCD64A1 burned area map (Fig. 9), also suggested that summer 2003 had 91.45 % more severe wildfire counts than summer 2010. It should be noted that high GFDI values do not guarantee a fire as there is no accounting for ignition sources, rather a higher GFDI value indicates higher probability of fire ignition and that if a grassland fire were to start it would spread faster compared to low GFDI values, given no fire suppression activity. Further complicating comparison of Fig. 8 and Fig. 9 is the presence of prescribed burns 10 that are deliberately done during low to moderate GFDI conditions, and that some fires shown in Fig. 9 occur in forested areas where GFDI is not applicable. Nevertheless, they provide a picture of the inter-annual spatial variability in both GFDI and burned area. The Weston Creek area is mostly comprised of forest with mixed land cover, whereas the Toodyay area is mostly a mix between croplands and savannas.
Using a burned area observation dataset from MODIS (MCD64A1), we test the effectiveness of GFDI with variable curing in increasing the probability that fires will occur in high GFDI severity levels compared to the probability that fires will occur in 25 low-moderate GFDI severity levels. Low intensity fires, such as prescribed burned, are removed from the burned area observation by using the FRP provided in MODIS active fire product (MCD14ML) to mask out burned area that have low FRP. We assume any burned area with FRP lower than 100 MW to be unburned (associated with low-moderate GFDI risk).
At each burned and unburned daily data point, the corresponding daily GFDI was calculated. The GFDI histogram in Fig. 11 shows the frequency of satellite based recalculated GFDIs and constant based (DOC = 100 %) reference GFDIs over burned 30 and unburned areas. Figure 11 shows that the recalculated GFDI places the largest percentage of unburned pixels in the lowmoderate GFDI severity class, with ~75 % of all unburned pixels occurring in the low-moderate or high severity classes.
Meanwhile the reference (DOC = 100 %) GFDI places ~80 % of unburned pixels in the high, very high, severe and extreme classes.
We can evaluate the performance in correctly assigning burned and unburned area for both recalculated and reference GFDI 35 by using the concept of ROC, as described earlier. Assume that the MCD64A1 burned area map represents the true condition and that the GFDI severity level represents the predicted condition, where the prediction is positive when GFDI level is classified as high or above for a burned area and low-moderate for an unburned area. Table 3 shows the contingency table, including both type I (unburned area with high or above GFDI level; false positive) and type II (burned area with low-moderate GFDI level) errors. Though recalculated GFDI has a lower true positive rate of correctly assigning burned area than reference 40 GFDI (0.86 vs 0.95), it is much better at assigning unburned area correctly, i.e. lower false positive rate (0.38 vs 0.53). Overall accuracy for recalculated GFDI is higher than the reference GFDI (0.62 vs 0.47).
Both Method B and MapVic DOC are then used to compute recalculated GFDI and compare with burned area observation dataset in the same manner as our DOC model. From the ROC analysis in Table 4 for Method B and MapVic recalculated GFDI, we found that even though Method B has the best DOC evaluation results (highest r 2 , lowest RMSE) and highest overall recalculated GFDI burned and unburned detection accuracy at 0.84, it is the worst at detecting burned area correctly with a true positive rate of only 0.10. This concurs with the findings in Newnham et al. (2010) who found Method B to consistently 5 under predict DOC, and hence it produces fewer cases of high or above GFDI severity. Our model is in the middle ground between Method B and MapVic in terms of overall accuracy.

Evaluation of Satellite based DOC
Previous studies that derived satellite based DOC have mostly relied solely on NDVI as a predictor. In the study by Newnham 10 et al. (2011) various forms of NDVI were used, including a straight NDVI linear regression and relative NDVI, as shown in Eq. (2). Their results suggested that a linear regression model based on NDVI alone reproduced site observations with an r 2 of 0.50. The results presented here indicate that inclusion of VOD in the regression model yields a comparable performance with r 2 of 0.44 to 0.55 (with independent sites evaluation at the worst and all sites evaluation at the best performance). However, when we compare the evaluation results with the currently available satellite based DOC from Bureau of Meteorology, the 15 best performer is Method B in both r 2 and RMSE for both evaluation methods (Newnham et al., 2010). We note that this is not a fair comparison as all the data used to evaluate the models was used to calibrate Method B, regardless of our evaluation methods. MapVic , on the other hand, was developed using a totally independent visual assessment dataset and hence is being evaluated against an independent dataset, regardless of our evaluation methods. Our first model performs better in all sites evaluation (include subset of the data for calibration) than independent sites evaluation, as expected. However, 20 our second model performs better in independent sites evaluation than all sites evaluation. Regardless, the overall performance is still not obviously better than neither Method B or MapVic models. While adding VOD to the DOC estimation model may introduce a comparable alternative to the existing optical based models, there is still no clear advantage of including VOD over the current methods.
An earlier study for estimating DOC directly with NDVI yielded even smaller RMSE of up to 6.3 %, but that particular study 25 is focused on data from only three different sites, within a limited study area of 1 km 2 (Dilley et al., 2004). Older studies that have used NDVI data derived from the National Oceanic and Atmospheric Administration's (NOAA) Advanced Very High Resolution Radiometer (AVHRR) have the same problem of interference from clouds and atmospheric effects, but do not have the advantage that VOD is not affected by clouds or aerosol interference.

GFDI with Various DOC estimates 30
Though the overall recalculated GFDI from Method B DOC is the best (overall accuracy from best to worst is 0.85 for Method B, 0.67 for our DOC model, 0.61 for MapVic, and 0.48 for GFDI with 100 % constant DOC), we found that it is the worst at detecting burned area correctly (true positive rate from best to worst is 0.88 for GFDI with 100 % constant DOC, 0.74 for MapVic, 0.69 for our DOC model, and 0.09 for Method B). Our VOD and NDVI based DOC model (first model) has a good balance in having the second best evaluation result and overall recalculated GFDI accuracy with a decent correct burned area 35 detection rate. We note that GFDI is only an indicator for the level of fire risk, and does not guarantee that a fire will occur, even at an extreme danger level. However, the improvement in accuracy indicates that inclusion of time and space varying DOC estimates makes it much more likely that areas identified at a GFDI severity level of high or above will burn than a lowmoderate severity level area.

Limitations
It is worth noting here that in an operational setting atmospheric interference by clouds or smoke will cause gaps in the optical (NDVI) data, though the VOD data remains unaffected. We also note that while the VOD data used here was derived from the AMSR-E sensor, which is no longer operational, VOD data derived from currently operating passive microwave sensors, such as the Advance Microwave Scanning Radiometer 2 (AMSR2), could be used in an operational setting. It should also be noted 5 that VOD's moderately coarse resolution of 0.1° may not be fine enough for use in many applications.
Reducing the chance of incorrectly assigning unburned and burned areas correctly from the ROC analysis made here is purely based on using the burned area map as a true baseline. However, the burned area map may include fires that are deliberately lit in low-moderate conditions, such as prescribed burns and fires that the GFDI is not designed for, such as a fire that burns in forested regions. Prescribed burns and low intensity fires are however minimised by applying low FRP threshold, using 10 information from the MODIS active fire product. The ROC analysis result here is only used to reinforce the idea that using the reference GFDI with constant curing (100 %) leads to overestimating GFDI in some situations, and might result in misleading fire danger warnings.
The satellite based DOC produced here is also at a moderate spatial resolution, which is a limitation of many satellite products.
However, DOC in reality can vary over spatial scales much finer than the satellite footprint (less than 500 m). As such, our 15 model should only be used as a guide for dynamic, near daily assessment of grassland curing at coarse to moderate spatial scales. This is also true for other satellite based DOC models, including Method B and MapVic models.

Conclusions
This study developed an alternative approach for estimating the grassland DOC using a relationship between the observed DOC and satellite based VOD and NDVI. The satellite based dataset was evaluated against the observed DOC data, which 20 resulted in a comparable performance with the currently existing optical based DOC estimation models. Despite the relatively coarse spatial resolution and temporal coverage of VOD and NDVI datasets used in this study, the satellite based DOC dataset produced from our model has the potential to contribute to the preparedness of fire management agencies and improve fire spread modelling. With a comparable, and arguably more balanced, performance in correctly predicting burned and unburned area through GFDI than currently available satellite based DOC models (i.e. Method B and MapVic), our model could provide 25 an appealing alternative estimated DOC data for GFDI computations and fire risk modelling.

Data availability
The following datasets and their associated sources or contact points are as listed below:

Author contribution
Wasin Chaivaranont, Jason Evans, and Yi Liu are involved in initial research planning and experiments design. For the rest of the study, including the analysis of the results, paper writing, and manuscript proofreading, all authors (Wasin Chaivaranont, 10 Jason Evans, Yi Liu, and Jason Sharples) are involved.

Competing interests
The authors declare that they have no conflict of interest.

Acknowledgements
We would like to acknowledge the following parties: 15 • NASA for making MODIS data freely available.
• NCI and CSIRO for making a freely available preprocessed MODIS data for Australia.
• Bureau of Meteorology for making freely available satellite based DOC dataset and continue to generate and distribute the data products set up by the AWAP project.
• Bushfire and Natural Hazards Cooperative Research Centre and its partner agencies for allowing us to use their 20 observed DOC dataset from their legacy project.
• AWAP for making daily Australia climate gridded dataset freely available.
• ECMWF for making ERA-Interim reanalysis gridded dataset freely available.