Measuring savanna woody cover at scale to inform ecosystem restoration

. The maturity of remote sensing and ecosystem restoration science provides new opportunities to monitor and assess ecosystem indicators at ﬁ ner resolutions and at suitable scales. A key link in joining these ﬁ elds is creating a framework to apply these rich datasets to ecosystem restoration projects. Savanna woodland ecosystems have discontinuous tree cover that is an important, but potentially highly variable, structural component of these ecosystems. Woody cover in savannas is also well suited to being measured with remote sensing techniques. We extracted woody cover from a 66-yr time series of aerial imagery covering 1718 ha of mesic savanna ecosystem in northern Australia, adjacent to the Ranger uranium mine and encompassing portions of Kakadu National Park. This ecosystem is proposed as a reference ecosystem that may be used to comparatively assess and monitor restoration trajectories of the mine site in the coming decades. The spatiotemporal patterns of change in woody cover were assessed at spatial extents similar to the mine site. We were able to construct a robust distribution of canopy cover values and associated spatial heterogeneity that can be used to set closure criteria, inform restoration trajectories, and guide monitoring activities for the restored mine site.


INTRODUCTION
A prevailing method for the assessment of ecosystem restoration projects requires the identification and characterization of a control or a desired end state (Wortley et al. 2013). In many cases, a proxy for the desired end state is used and is often called a reference ecosystem. A reference ecosystem will contain compositional, structural, and functional attributes that are measurable (McDonald et al. 2016), and at some spatiotemporal scale, predictable. These attributes provide the basis for the development of specific indicators which are used to set goals for restoration projects and to assess progress toward recovery and achievement of the desired restored ecosystem (McDonald et al. 2016).
In many restoration projects, including the restoration of mine sites, quantitative indicators representing attribute components (e.g., nutrient cycling, vegetation cover, species diversity) are often used as closure criteria, which are a set of conditions that formally define restoration success. Indicators form the basis of restoration monitoring programs and may inform management strategies (Thom 2000, Convertino et al. 2013. The selection of indicators requires important consideration. Aside from the representativeness of an indicator for a specific attribute, the spatiotemporal variability of an indicator should be understood before it is developed and used in a restoration framework, especially (1) the relationship between the spatial scale of indicator variability and the physical size of a restoration project, and (2) the temporal variability of an indicator and the time horizon of a restoration project. It may be difficult to partition indicator variability into process-driven (e.g., spatially explicit fire history in a savanna) and stochastic (e.g., windthrow from storm events) categories as large-scale experiments designed for this purpose may be impractical. This should not dissuade us from identifying what underlying ecological processes an indicator may represent, and how variability or patchiness in that indicator informs our understanding of the structure and function of an ecosystem (Ludwig and Tongway 1995). This leads to an important question: "At what spatial and temporal scales should we measure an indicator across reference ecosystems and restoration sites?".
For planning and design purposes, restoration practitioners require knowledge of the likely range of values of an indicator before the project even commences, and well in advance of the time that closure criteria are expected to be met. In some instances, there is practical value in developing long time series of indicators from reference ecosystems in order to better define the range in values of the indicators and any underlying trends (e.g., climate change influences; Large et al. 2007).
Covering over one-eighth of the Earth's land area, savannas are one of the largest terrestrial biomes (Scholes andArcher 1997, Sankaran et al. 2004). The two main structural components of savanna are a canopy of discontinuous tree cover above a mostly continuous grass understorey (Hutley and Setterfield 2008). The tree-grass ratio in savanna varies over time and space (Stevens et al. 2017), but there remain significant knowledge gaps regarding the drivers of tree-grass dynamics in these ecosystems (Staver et al. 2017). At a landscape scale, the tree-grass dynamics may be driven by soil-type, nutrient availability, fire, herbivory, extreme weather events (such as cyclones), and human activity (Staver et al. 2017). In regions where the mean annual precipitation is above 650 mm, fire and herbivory may play an important role in maintaining savannas and preventing woody canopy closure (i.e., forming forest; Sankaran et al. 2005. With a range of drivers and their interactions, it is not surprising that any predictive model of heterogeneity in these ecosystems remains poorly resolved. Staver (2018) outlines how at larger spatial extents (potentially those relevant to large-scale ecosystem restoration), tree density and demographics are the most variable components of the ecosystem, whereas at these scales, grass cover and biomass are reasonably predictable. They further explain how sufficient sampling and choice of plot sizes should be carefully considered in order to adequately describe the heterogeneity in savanna tree layers.
When considering indicators of restoration success that can be measured over long time periods, it is important to consider the context of local, regional, and global trends that may be influencing these ecosystems. On a regional scale pertinent to this study, there are several longterm studies of vegetation dynamics in Australian tropical savanna; however, these tend to be limited in spatial context-either using noncontiguous data or sub-samples of an area , or limited in temporal resolution-using only two dates (Bowman et al. 2001, Sharp andBowman 2004). The lack of spatial and temporal contiguity can be attributed to data availability and limited computer processing resources at the time of study. Lehmann et al. (2008) found savanna tree cover varied greatly for their sites in the Kapalga region in northern Australia between 1964 and 2004. The dynamism of tree cover was driven by the frequency and intensity of fire, which in turn can be influenced by the grazing of feral animals, specifically water buffalo, which suppressed perennial species of grasses and the recruitment of trees to the canopy (Werner et al. 2006). Tropical cyclones and other extreme weather events are also potential agents of change (Lehmann et al. 2009).
In savannas of northern Australia, there has been reported thickening of woody cover (Chen et al. 2003, Banfai and Bowman 2005, Lehmann et al. 2009), while other studies report a reduction in savanna woody cover (Sharp and Bowman 2004). Most of these studies are limited in either spatial coverage or temporal resolution. Recent research shows that there may be no thickening of savanna woody cover, with fire frequency and intensity preventing expansion through suppression of woody growth (Murphy et al. 2014, Prior et al. 2020. Remote sensing data, including historical air photograph archives and satellite imagery, can provide information for the analysis of landscapes over time and potentially at a resolution sufficient to derive appropriate indicators for large-scale ecosystem restoration. There is increasing focus on developing indicators that can be assessed at the scale of an entire restoration project (Doren et al. 2009). One indicator that can be successfully extracted from such data is woody cover (defined here as the canopy cover of the tree layer). Scale can play a role in describing variability of woody cover (Staver 2018). Metrics for whole of site or large patches can poorly reflect or misrepresent patterns or processes operating at a small scale (Zhang et al. 2019). As such, measuring woody cover at the whole of site level may not capture the variability inherent in the landscape. Spatial heterogeneity is best represented at a scale where it is readily detected and measured.
The Australian Government's Supervising Scientist has developed an ecosystem restoration standard for the Ranger uranium mine in the Northern Territory of Australia to be used to (1) define indicators and associated metrics for similarity and sustainability attributes against which to assess success of restoration at the site, (2) monitor progress of ecosystem restoration against a restoration trajectory, and (3) assess achievement of closure criteria . The ecosystem restoration standard addresses the regulatory environmental requirement to establish an ecosystem that is similar to adjacent areas of Kakadu National Park (Department of the Environment and Heritage 1999). Woody cover is one of the indicators of the structural attributes for similarity from the standard. Although it would be preferable to compare reference ecosystem woody cover with that of the restored mine site, restoration activities, including revegetation, are just commencing and will continue for at least the next decade. As such, we have mapped and quantified woody cover over 10 dates from a proposed savanna reference ecosystem in the region surrounding Ranger to provide information on the dynamics of savanna woody cover across space and time. The analysis undertaken addresses the following questions: (1) How variable is woody cover over multiple spatial and temporal scales? (2) How do we incorporate variability of indicators in ecosystem assessments of large mine sites?
The answers to these questions will provide a foundation to enable an evidence-based comparison of reference and restored ecosystems across temporal and spatial scales.

Study area
The study area is described in detail in Whiteside et al. (2020) and is defined by a 10-km buffer around the disturbed footprint of the Ranger uranium mine, including parts of the UNESCO World Heritage-listed Kakadu National Park (Fig. 1). The vegetation of the region is predominantly a savanna matrix of eucalypt woodland with annual grass understorey, intersected by drainage lines dominated by Melaleuca spp. Soils in the area are mostly sandy loams formed over the underlying, weathered sandstone Koolpinyah surface (Wells 1979). Wells (1979) described the landforms in the region and categorized a range of slopes, soil types, and overlying vegetation into distinct land units. Four of the open forest/savanna woodland (non-riparian) Wells land units were selected for this study on the basis of the presence on the pre-mine site or location directly adjacent to the mine site as well as having slopes similar to the proposed restored mine site ( Table 1).
The region has a tropical monsoonal climate with the majority (over 90%) of rainfall falling in the months of November to April. Mean rainfall for nearby Jabiru Airport is 1565 mm (Bureau of Meteorology 2019). The mean maximum temperature is 34.4°C, while mean minimum is 22.7°C. Most of the annual grasses begin to cure at the beginning of the dry season creating a fuel load. Fires are common in the dry season (May-October; Edwards et al. 2003), and current management subscribes to small-scale patch burning in the early dry season to minimize the risk and impact of large and severe fires later in the dry season (Director of National Parks 2016). The western side of the study area has experienced more frequent fires than the eastern side .

Image data and processing
Image data covering the study area were available for 10 dates (Appendix S1: Table S1). The first eight dates of image data were mosaics of aerial photography, ranging from gray scale (1950 and 1964), true color (1978, 1981, 1984, 1987, and 2004), and near-infrared false color (1976). The last two datasets, 2010 and 2016, were high spatial resolution WorldView-2 and WorldView-3 satellite imagery, respectively. The spatial extent of each image set is shown in Appendix S1: Fig. S1. The high level of temporal  resolution over an extended period (10 epochs spanning 66 yr) is exceptional, with most published studies analyzing change from historical imagery limited to as few as two dates (e.g., Sharp and Bowman 2004, Platt and Schoennagel 2009, Levick and Rogers 2011. Image processing to extract woody cover was based upon Whiteside et al. (2020) but was expanded from 4 to 10 dates of image capture. The technique is best applied to data where the phenological difference between trees and understorey is maximized, that is, dry season capture, when understorey grasses have cured and tree canopies are still green.

Selection of proposed reference ecosystem
The 4 Wells land units have similar slopes to those proposed for the final constructed landform of the mine site and contain savanna woodland vegetation communities that dominate the adjacent areas of Kakadu National Park. Once image processing was complete and woody cover extracted, areas not contained in the 4 land units, roads, the disturbed areas of the mine site, and the town of Jabiru were removed from the dataset using ESRI's ArcGIS Pro. This created 59 separate polygons of the 4 different land units containing a total of 1718 hectare cells of woody cover data (Fig. 1). These polygons contain our basis for the reference ecosystem suitable for deriving an array of indicator metrics for assessing ecosystem restoration success at Ranger uranium mine.

Analysis
Analyses examined change in the proportion of woody cover at various spatial scales over time. The scales were per land unit type (n = 4), land unit polygon (n = 59), and by dividing the data into spatial units similar to the size (761 ha, n = 25) of the disturbed footprint of the mine site (see spatial subdivision description below). The spatial coverage of the 1978 aerial mosaic and the 2010 satellite imagery was not as complete as for the other dates, so areas where there were null values for either of those two years were removed from the time series analysis. We examined the per-polygon data in order to evaluate the effects of geography and time on woody cover using a linear model. The UTM coordinates of the centroid of each land unit polygon were calculated using ESRI's ArcGIS Pro software and were used to represent the geographic position of each polygon. To assess whether there was a geographically determined effect on woody cover, each land unit polygon and its coordinates were fit with a linear model with year specified as a factor.
The per hectare woody cover data for each land unit were aggregated over time, and a series of 2-sample Kolmogorov-Smirnov tests were conducted on all combinations of the 4 land units to test for differences in the distributions between land units. The same data were also analyzed using a one-way ANOVA.
The study area was divided into 25 approximately equal-sized blocks. The maximum latitudinal (13.4 km) and longitudinal (14.2 km) extents were divided by 5 and used to create the dimensions of each block (2.68 9 2.84 km, i.e., 761 ha; Appendix S1: Fig. S2). These blocks (761 ha) closely approximate the area requiring restoration at Ranger uranium mine. For each of the 25 blocks of savanna, the amount of woody cover was calculated from the relevant land units (with disturbed areas masked out) as previously described. This means that for each block of savanna woodland, only a portion of it was measured for woody cover. It is likely that ecosystem restoration at Ranger will primarily use a reference ecosystem of the vegetation community represented by the land units used in this study. The total area of the Ranger site with that relevant vegetation community will likely be larger than the total area represented in any of the 25 land blocks. Some blocks contained very small amounts of relevant vegetation and woody cover data. In order to minimize any effects from small sample sizes, blocks with <40 ha of data were excluded because they demonstrated withinblock variance characteristics not reflective of variances of blocks with more data available. This left 15 blocks containing between 42 and 242 ha of woody cover data from relevant land units. The means and within-block standard deviations of these blocks were assembled into a distribution. In order to understand the shape of these distributions, they were tested for normality using Shapiro-Wilk's test.
All statistical analyses were performed using the R programming language (R Core Team 2019) with figures produced using the ggplot2 package (Wickham 2016).

RESULTS
Since 1950, there has been no apparent change in the proportion of woody cover across the study site (Fig. 2). The minimum mean woody cover in a year is 27.6% per ha cell (1978), and the maximum is 43.0% (2010). The overall mean woody cover for all ha cells over time is 35.1% (AE 15.1% SD).
Across the study site, there is a slight geographic effect of increase in woody cover from the south to the north direction (Appendix S1: Fig. S3; Northing: P = 0.06, Adjusted r 2 = 0.017) and decrease from the east to west direction (Appendix S1: Fig. S4; Easting: P < 0.0001, Adjusted r 2 = 0.05). Thus, a negligible amount of the variance in woody cover is explained by geographic location across the study site.
After aggregating all time point data by the type of underlying land unit, a series of pairwise 2-sample Kolmogorov-Smirnov tests found differences in the distribution of per hectare woody cover values among all four land units (largest P = 0.0006; Table 2). When aggregating all time points, the minimum mean per ha proportion cover was 31.6% for land unit 4a, while the maximum was 39.4% for land unit 4c1 ( Table 2). The range in standard deviation between land units is minimal (0.140-0.164). Examining each land unit by time point (Fig. 3)

revealed inconsistent patterns between land units across years.
When dividing the study area into 25 spatial blocks of similar size to the proposed Ranger restoration site, there is a relationship between the standard deviation in woody cover and the number of 1-ha plots represented within that spatial block (Fig. 4), especially at small sample sizes. With larger sample sizes within a spatial block, the standard deviation begins to reach an asymptote close to the standard deviation of the entire dataset, although even with the largest sample size (n = 242 ha), the asymptote is not reached. Because of the potential for small sample sizes to influence the variability found within a spatial block, the next set of analyses was restricted to those spatial blocks with at least 40 ha of data (n = 15 blocks). The mean percentage of woody cover for these blocks is 35.8% across all time points (range 14.0-59.4%). The distribution of woody cover from these spatial blocks conforms to a normal distribution (Shapiro-Wilk's test, P = 0.176), with quartile values of 29.9% and 40.7%. When color coding these spatial blocks by year (Fig. 5), patterns in the data can be visualized. For example, the 2016 data have a narrower range of woody cover There is a significant but weak positive relationship between the mean canopy cover and standard deviation within a land block (Fig. 6, Pearson R = 0.27, P = 0.0007), but no apparent relationship between the area of reference ecosystem sampled and either the mean or standard deviation (Fig. 6). The per-spatial block standard deviations in canopy cover are not distributed normally (Shapiro-Wilk's test, P < 0.001) and are right skewed (Appendix S1: Fig. S5). There are similar interannual patterns in standard deviation to those found in the distribution of the means of woody cover, with the 1950-2016 data comparison similarly showing a narrower range of woody cover standard deviation in 2016 than 1950.

Variability of woody cover over spatial and temporal scales
The context for addressing our first research question, "How variable is woody cover over multiple spatial and temporal scales?", is that savanna ecosystems have spatiotemporal patterns and drivers and stressors that can make these systems patchy and dynamic. As a result, the scale of measurement may strongly influence the interpretation of how savannas are described (Price et al. 2005, Russell-Smith and Yates 2007, Riginos et al. 2009, Staver 2018. The spatial resolution and temporal extent required to fully characterize these   ecosystems and their responses to environmental drivers need to be carefully considered (Mayer and Rietkerk 2004, Staver et al. 2017, Staver 2018. Additionally, the stressors on, and drivers of, savanna ecosystems operate across a large range of scales to which woody cover may behave nonlinearly through space and time (Sankaran et al. 2005, Wigley et al. 2010, Hutley et al. 2011, Murphy and Bowman 2012, Staver et al. 2017. Measuring indicators of ecosystem restoration at spatial extents similar to the size of a restoration project (referred to as at-scale) and over decadal time periods allows for a greater understanding of their variability. This approach improves confidence about how an indicator may perform when used in a quantitative assessment.
There are slight spatial trends of woody cover in the cardinal directions, with the east to west direction having a weak but consistent trend through time that may be related to fire history and frequency . Ranger is positioned close to the center of the surrounding reference ecosystem, and as a result, small geographic differences may be minimized by developing indicator values from all available directions drawn from reference ecosystem within 10 km of the mine site. Additionally, the land units identified as containing the reference ecosystem have similar woody cover. It is important to note that this study does not address the drivers of variation in woody cover, as Prior et al. (2020) concluded that fire history, recent rainfall, and image capture date did not have any influence across the study site over the same time span.
By dividing the reference ecosystem into spatial blocks approximately the same size as the mine site, we can aggregate the data in a way that may be more representative of what could be expected for a site such as Ranger. An important consideration in this process is to evaluate the relationship between sample size and patchiness within a spatial block. In this way, when using at-scale metrics of indicators as an assessment tool, it is appropriate to expect that the metrics for the respective restoration site indicator fit in the envelope of variability defined for the broader reference ecosystem. There may be a threshold where spatial blocks do not contain adequate reference ecosystem to be used as an indicator as a result of failing to capture the heterogeneity measured in spatial blocks containing more reference ecosystem.

Application to ecosystem restoration of large mine sites
With regard to the research question "How do we incorporate variability of indicators in ecosystem assessments of large mine sites?", the inherent variability of canopy cover needs careful consideration. In some regards, the local scale patchiness in reference ecosystems may allow for a wide envelope of what is considered successful restoration. On the other hand, it appears that patchiness is itself a property of savanna canopy cover. By measuring an indicator at the same scale as a restoration project and developing appropriate replicates of the measured indicator, a within-and between-spatial block variance can be calculated. The distribution of within-block variance can be used as an estimate of woody cover heterogeneity that is observed within spatial blocks approaching the size of a restoration site. It is not surprising that the within-block variance in woody cover is lower than the overall variance of woody cover.
We propose pairing the approach of spatial blocking indicator data with a risk analysis framework ensuring that the variability of the broader ecosystem is not used to unduly influence expectations and improperly inform mine closure criteria. Thus, while the tails of the distribution (Fig. 5) fall within the bounds of woody cover that is seen in comparatively sized reference ecosystem spatial blocks, they may also represent risk of failed restoration (e.g., trajectory deviation, potentially interpreted as over-or under-planting). Measuring indicators at appropriate spatial and temporal scales improves parameterization of the patterns of central tendency and heterogeneity of a reference ecosystem and best serves a risk-minimization approach.

Trajectories and monitoring
The drivers and stressors in savannas provide important context for understanding the trajectory of a restoration project in this type of ecosystem. Furthermore, monitoring ecosystem development in a trajectory context needs well-defined indicators that can be readily measured and used to inform management processes (Ludwig et al. 2003, Tongway andLudwig 2006). Through effective monitoring, in this case the regular capture of aerial or satellite imagery, a restoration site can be compared through time with a paired reference ecosystem. At spatial extents similar to the rehabilitation of Ranger mine, a proposed savanna woodland reference ecosystem has canopy cover with a large range when measured at a plot size of 1 ha.
When randomly subsampling the data, a smaller number of 1-ha plots is often able to characterize the entire dataset of woody cover at each time point. In a monitoring program, this may be useful information if data collection methods are costly or time consuming. However, it is important to note that this practice would not capture the spatial heterogeneity or patchiness expected in a restoration site, as the subsampling would need to be conducted on spatial blocks of similar size to the restoration site; hence, if a goal of restoration is to incorporate the spatial heterogeneity measured in reference ecosystems, this approach may not be advisable.
Canopy cover can be useful as a structural indicator to monitor how rehabilitation is progressing over time because it can be measured using cost-effective, repeatable techniques across both an entire restoration site and its reference ecosystem. Furthermore, spatial blocks can be delineated to further assess the variability that is occurring across a reference ecosystem and used to assess the status of a restoration site. As a result, if the restored ecosystem is on a desired trajectory, canopy cover values should approach the range of at-scale values as measured over a reference ecosystem.

Assessment methods and closure criteria
Conventional methods for assessing mine site rehabilitation often deal with threshold (canopy cover per ha must be between X% and Y%) or means-based comparative (canopy cover per hectare must have a mean similar to that found in the reference ecosystem) approaches. Assessment of canopy cover, like a number of other heterogeneous structural attributes of savanna ecosystems, may not be amenable to these approaches, or at least is potentially prone to unforeseen outcomes if not prescribed appropriately in closure criteria. Threshold-and means-based approaches tend to neglect the inherent variability of the indicator. As a result of the spatial variability demonstrated in canopy cover, an unforeseen outcome of means or threshold approaches may be a restored site that is more homogenous than the reference ecosystem. Comparing measures of variability of canopy cover values between ecosystems may best imbue the heterogeneity observed and expected within a restored site when compared to a reference ecosystem. When data can be captured at large spatial extents, a two-sample Kolmogorov-Smirnov test and its derivatives may be candidate statistical tests for comparing two distributions (a reference and restored ecosystem).
If an indicator can be replicated at spatial extents that are similar to the restoration site, a distribution of central tendency and a distribution of the variance (representing heterogeneity) of spatial units from within the reference ecosystem can be constructed. In this case, a means or threshold approach can be applied to both the central tendency and variation of the indicator and applied in a closure framework. In this scenario, the closure framework should employ a parallel risk analysis that may conclude that the tails of the distribution represent ecosystems that may not be on a suitable trajectory and have not met closure criteria. Identifying what portion of the tails of the distribution of sparser and thicker woody cover, respectively, do not meet closure criteria may be determined by analyzing the range of variability in woody cover that individual spatial blocks undergo through time. If cover in reference ecosystem spatial blocks commonly wax and wane between sparse cover to moderate cover or from thick woody cover to moderate cover over time, there may be more confidence in including data from the entire distribution, that is, including data from the tails of the distribution, as having met closure criteria because the behavior of individual spatial blocks of reference ecosystem woody cover regularly approaches the measure of central tendency from the entire reference ecosystem.

Summary
As technology increases the extent and resolution of data capture, our capacity to monitor and measure ecosystems at suitable scales improves. As scalability improves, we need to carefully examine the variability and heterogeneity that is measured from ecosystem attributes with regard to how they are applied to the science of ecosystem restoration. It is important to understand whether and how these technological improvements can be leveraged toward better restoration outcomes. By measuring woody cover at-scale relevant to large mine site ecosystem restoration projects, the data can be used to inform restoration practices throughout the life cycle of a project. Importantly, at-scale indicators such as woody cover can be used to fully describe the heterogeneity seen through space and time in comparable reference ecosystems.